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CN Abstract 

^~i In this paper we investigate a potential use of fluid approximation techniques in the context 

V^ of stochastic model checking of CSL formulae. We focus on properties describing the behaviour 

)—^ of a single agent in a (large) population of agents, exploiting a limit result known also as 

(/3 fast simulation. In particular, we will approximate the behaviour of a single agent with a 

^ O^ time-inhomogeneous CTMC, which depends on the environment and on the other agents only 

through the solution of the fluid differential equation, and model check this process. We will 

CNJ prove the asymptotic correctness of our approach in terms of satisfiability of CSL formulae. We 

K*" will also present a procedure to model check timc-inhomogcneous CTMC against CSL formulae. 

(*-5^ Keywords: Stochastic model checking; fluid approximation; mean field approximation; reach- 

(^ ability probability; time-inhomogeneous Continuous Time Markov Chains 
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In recent years, there has been a growing interest in fluid approximation techniques in the formal 
methods community [U El EJ SJ El |6] . These techniques, also known as mean field approximation, 
are useful for analysing quantitative models of systems based on continuous time Markov Chains 
(CTMC), possibly described in process algebraic terms. They work by approximating the discrete 
state space of the CTMC by a continuous one, and by approximating the stochastic dynamics of 
the process with a deterministic one, expressed by means of a set of differential equations. The 
asymptotic correctness of this approach is guaranteed by limit theorems [3 El E]) showing the 
convergence of the CTMC to the fiuid ODE for systems of increasing sizes. 

The notion of size can be different from domain to domain, yet in models of interacting agents, 
usually considered in computer science, the size has the standard meaning of population number. 
All these fluid approaches, in particular, require a shift from an agent-based description to a 
population-based one, in which the system is described by variables counting the number of agents 
in each possible state and so individual behaviours are abstracted. In fact, in large systems, the 
individual choices of single agents have a small impact, hence the whole system tends to evolve 
according to the average behaviour of agents. Therefore, the deterministic description of the 
fluid approximation is mainly related to the average behaviour of the model, and information 



about statistical properties is generally lost, although it can be partially recovered by introducing 
fluid equations of higher order moments of the stochastic process (moment closure techniques 

uniiiiKis]). 

Differently to fluid approximation, the analysis of quantitative systems like those described 
by process algebras can be carried out using quantitative model checking. These techniques have 
a long tradition in computer science and are powerful ways of querying a model and extracting 
information about its behaviour. As far as stochastic model checking is considered, there are some 
consolidated approaches based mainly on checking Continuous Stochastic Logic (CSL) formulae 
[I3l mi [15], which led to widespread software tools [16]. All these methods, however, suffer (in a 
more or less relevant way) from the curse of state space explosion, which severely hampers their 
practical applicability. In order to mitigate these combinatorial barriers, many techniques have 
been developed, many of them based on some notion of abstraction or approximation of the original 
process pTlfTS]. 

In this paper, we will precisely target this problem, trying to see to what extent fluid approx- 
imation techniques can be used to speed up the model checking of CTMC. We will not tackle 
this problem in general, but rather we will focus on a restricted subset of system properties: We 
will consider population models, in which many agents interact, and then focus on the behaviour 
of single agents. In fact, even if large systems behave almost deterministically, the evolution of 
a single agent in a large population is always stochastic. Single agent properties are interesting 
in many application domains. For instance, in performance models of computer networks, like 
client-server interaction, one is often interested in the behaviour and in quality-of-service metrics 
of a single client (or a single server) , such as the waiting time of the client or the probability of a 
time-out. 

Single agent properties may also be interesting in other contexts. For instance, in ecological 
models, one may be interested in the chances of survival or reproduction of an animal, or in its 
foraging patterns [19| . In biochemistry, there is some interest in the stochastic properties of single 
molecules in a mixture (single molecule enzyme kinetics [201 IZI] ) • Other examples may include the 
time to reach a certain location in a traffic model of a city, or the chances to escape successfully 
from a building in case of emergency egress [22] . 

The use of fluid approximation in this restricted context is made possible by a corollary of 
the fluid convergence theorems, known by the name of fast simulation [23', ^ , which provides a 
characterization of the behaviour of a single agent in terms of the solution of the fluid equation: 
the agent senses the rest of the population only through its "average" evolution, as given by the 
fluid equation. This characterization can be proved to be asymptotically correct. 

Our idea is simply to use the CTMC for a single agent obtained from the fluid approximation 
instead of the full model with A^ interacting agents. In fact, extracting metrics from the descrip- 
tion of the global system can be extremely expensive from a computational point of view. Fast 
simulation, instead, allows us to abstract the system and study the evolution of a single agent 
(or of a subset of agents) by decoupling its evolution from the evolution of its environment. This 
has the effect of drastically reducing the dimensionality of the state space by several orders of 
magnitude. 

Of course, in applying the mean field limit, we are introducing an error which is difficult to 
control (there are error bounds but they depend on the final time and they are very loose [U]). 
However, this error in practice will not be too large, especially for systems with a large pool 
of agents. We stress that these are the precise cases in which current tools suffer severely from 
state space explosion, and that can mostly benefit from a fluid approximation. However, we will 
see in the following that in many cases the quality of the approximation is good also for small 



populations. 

In the rest of the paper, we will basically focus on how to analyse single agent properties of 
three kinds: 

• Next-state probabilities, i.e. the probability of jumping into a specific set of states, at a 
specific time. 

• Reachability properties, i.e. the probability of reaching a set of states G, while avoiding 
unsafe states U. 

• Branching temporal logic properties, i.e. verifying CSL formulae. 

A central feature of the abstraction based on fluid approximation is that the limit of the model 
of a single agent has rates depending on time, via the solution of the fluid ODE. Hence, the limit 
models are time-inhomogeneous CTMC (ICTMC). This introduces some additional complexity in 
the approach, as model checking of ICTMC is far more difficult than the homogeneous-time case. 
To the best of the author's knowledge, in fact, there is no known algorithm to solve this problem 
in general, although related work is presented in Section [2j We will discuss a general method 
in Sections |4j [5] and [6j based on the solution of variants of the Kolmogorov equations, which is 
expected to work for small state spaces and controlled dynamics of the fluid approximation. The 
main problem with ICTMC model checking is that the truth of a formula can depend on the time 
at which the formula is evaluated. Hence, we need to impose some regularity on the dependency of 
rates on time to control the complexity of time-dependent truth. We will see that the requirement, 
piecewise analyticity of rate functions, is intimately connected not only with the decidability of 
the model checking for ICTMC, but also with the lifting of convergence results from CTMC to 



truth values of CSL formulae (Theorems 6.1 and |6.3[ ). 

The paper is organized as follows: in Section [2] we discuss work related to our approach. In 
Section [3| we introduce preliminary notions, fixing the class of models considered (Sectio n |3.1[ ) 



and presenting fluid limit and fast simulation theorems (Sections |3.2 and 3.3). In Section I4| we 
present the algorithms to compute next-state probability. In Section 5j instead, we consider the 
reachability problem, presenting a method to solve it for ICTMC. In both cases, we also discuss the 
convergence of next-state and reachability probabilities for increasing population sizes. In Section 
[6j instead, we focus on the CSL model checking problem for ICTMC, exploiting the routines for 
reachability developed before. We also consider the convergence of truth values for formulae about 
single agent properties. Finally, in Section[7| we discuss open issues and future work. All the proofs 
of propositions, lemmas, and theorems of the paper are presented in Appendix [A| A preliminary 
version of this work has appeared in pi] . 

2 Related work 

Model checking (time homogeneous) Continuous Time Markov Chains (CTMC) against Contin- 
uous Stochastic Logics (CSL) specifications has a long tradition in computer science [T3 l [T H [T5]. 
At the core of our approach to study time-bounded properties there are similarities to that de- 
veloped in [13j , because we consider a transient analysis of a Markov chain whose structure has 
been modified to reflect the formula under consideration. But the technical details of the transient 
analysis, and even the structural modification, differ to refiect the time-inhomogeneous nature of 
the process we are studying. 

In contrast, the case of time-inhomogeneous CTMCs has received much less attention. To the 
best of the authors' knowledge, there has been no previous proposal of an algorithm to model 



check CSL formulae on a ICTMC. Nevertheless model checking of ICTMCs has been considered 
with respect to other logics. Specifically, previous work includes model checking of HML and LTL 
logics on ICTMC. 

In |25j . Katoen and Mereacre propose a model checking algorithm for Hennessy-Milner Logic 
on ICTMC. Their work is based on the assumption of piecewise constant rates (with a finite 
number of pieces) within the ICTMC. The model checking algorithm is based on the computation 
of integrals and the solution of algebraic equations with exponentials (for which a bound on the 
number of zeros can be found). 

LTL model checking for ICTMC, instead, has been proposed by Chen et al. in |26j . The 
approach works for time-unbounded formulae by constructing the product of the CTMC with a 
generalized Biichi automaton constructed from the LTL formula, and then reducing the model 
checking problem to computation of reachability of bottom strongly connected components in 
this larger (pseudo)-CTMC. The authors also propose an algorithm for solving time bounded 
reachability similar to the one considered in this paper (for time-constant sets). 

Another approach related to the work we present is the verification of CTMC against deter- 
ministic time automata (DTA) specifications [27J, in which the verification works by taking the 
product of the CTMC with the DTA, which is then converted into a Piecewise Deterministic 
Markov Process (PDMP, [28]), and then solving a reachability problem for the so obtained PDMP. 
This extends earlier work by Baier et al. [29J and Donatelli et al. [30j . These approaches were lim- 
ited to considering only a single clock. This means that they are albe to avoid the consideration of 
ICTMC, in the case of |30| . through the use of supplementary variables and subordinate CTMCs. 

In [31] , Chen et al. consider the verification of time-homogenenous CTMC against formulae in 
the the metric temporal logic (MTL). This entails finding the probability of a set timed paths that 
satisfy the formula over a fixed, bounded time interval. The approach taken is one of approxima- 
tion, based on an estimate of the maximal number of discrete jumps that will be needed in the 
CTMC, N, and timed constraints over the residence time within states of a path with up to N 
steps. The probabilities are then determined by a multidimensional integral. 

Our work is underpinned by the notion of fast simulation, which has previously been applied 
in a number of different contexts [9]. One recent case is a study of policies to balance the load 
between servers in large-scale clusters of heterogeneous processors [23] . A similar approach is 
adopted in [321, in the context of Markov games. These ideas also underlie the work of Hayden et 
al. in |33j . Here the authors extend the consideration of transient characteristics as captured by the 
fluid approximation, to approximation of first passage times, in the context of models generated 
from the stochastic process algebra PEPA. Their approach for passage times related to individual 
components is closely related to the fast simulation result and the work presented in this paper. 

3 Preliminaries 

In this section, we will introduce some backgound material needed in the rest of the paper. First 
of all, we introduce a suitable notation to describe the population models we are interested in. 



This is done in Section 3.1 In particular, models will depend parametrically on the (initial) 
population size, so that we are in fact defining a sequence of models. Then, in Section |3.2[ we 
present the classic fluid limit theorem, which proves convergence of a sequence of stochastic models 



to the solution of a differential equation. In Section 3.3 instead, we describe fast simulation, a 
consequence of the fluid limit theorem which connects the system view of the fluid limit to the 
single agent view, providing a description of single agent behaviour in the limit. Finally, in Section 



3.4, we recall the basics of Continuous Stochastic Logic (CSL) model checking. 



3.1 Modelling Language 

In the following, we will describe a basic language for CTMC, in order to fix the notation. We 
have in mind population models, where a population of agents, possibly of different kinds, interact 
together through a finite set of possible actions. To avoid a notational overhead, we assume that 
the number of agents is constant during the simulation, and equal to N. Furthermore, we do not 
explicitly distinguish between different classes of agents in the notation. 

In particular, let Y^ ^ S represent the state of agent i, where S = {1,2, ... ,n} is the 
state space of each agent. Multiple classes of agents can be represented in this way by suitably 
partitioning S into subsets, and allowing state changes only within a single class. Notice that we 
made explicit the dependence on A^, the total population size. 

A configuration of a system is thus represented by the tuple {Y^ , ■ ■ ■ , Y}^ ) . When dealing with 
population models, it is customary to assume that single agents in the same internal state cannot 
be distinguished, hence we can move from the agent representation to the system representation 
by introducing variables counting how many agents are in each state. With this objective, define 

N 

xf^=Y^l{Yl-^=j}, (1) 

i=l 

so that the system can be represented by the vector X'^^ = (X| , . . . ,Xn ), whose dimension 
is independent of N. The domain of each variable X- is obviously {1, . . . , A^}. 

We will describe the evolution of the system by a set of transition rules at this global level. 
This simplifies the description of synchronous interactions between agents. The evolution from the 
perspective of a single agent will be reconstructed from the system level dynamics. In particular, 
we assume that X^^' is a CTMC (Continuous-Time Markov Chain), with a dynamics described 
by a fixed number of transitions, collected in the set T^ ' . Each transition r S ']'("> is defined 
by a multi-set of update rules Rr and by a rate function r^ . The multi-seljj -Rr contains update 
rules p £ Rr oi the form i ^ j, where i,j £ S. Each rule specifies that an agent changes state 
from i to j. Let rrir^i^j denote the multiplicity of the rule i — )• j in R-^. We assume that R^ is 
independent of A^, so that each transition involves a finite and fixed number of individuals. Given 
a multi-set of update rules Rr, we can define the update vector v,- in the following way: 

where 1, is the vector equal to one in position i and zero elsewhere. Hence, each transition changes 
the state from X^ -* to X*^ •* + v,-. The rate function rj- (X) depends on the current state of the 
system, and specifies the speed of the corresponding transition. It is assumed to be equal to zero 
if there are not enough agents available to perform a r transition. Furthermore, it is required to 
be Lipschitz continuous. We indicate such a model by X^ > = (X^ ',?"' ^Xq ), where Xq is 
the initial state of the model. 

Given a model X^ ', it is straightforward to construct the CTMC associated with it, ex- 
hibiting its infinitesimal generator matrix. First, its state space is T> = {{xi, . . . ,Xn) \ Xi G 
{1, . . . ,N},^^Xi = N}. The infinitesimal generator matrix Q, instead, is the D x D matrix 
defined by 

gx,x' = '^{rr{:x.) I r G T, x = x + v^}. 

We will indicate the state of such a CTMC at time t by X(t). 



'^The fact that Rr is a multi-set, allows us to model events in which agents in the same state synchronise. 




Figure 1: Visual representation of the client server system of the running example. 



Example. We introduce now the main running example of the paper: we will consider a model 
of a simple client-server system, in which a pool of clients submits queries to a group of servers, 
waiting for a reply. In particular, the client asks for information from a server and waits for it 
to reply. It can time-out if too much time passes. The server, instead, after receiving a request 
does some processing and then returns the answer. It can time-out while processing and while it 
is ready to reply. After an action, it always logs data. The client and server agents are visually 
depicted in Figure [T} The global system is described by the following 8 variables: 

• 4 variables for the client states: Crq, C^, Cj-c, and Ct- 

• 4 variables for the server states: Sj-q, Sp, Srp, and Si. 

Furthermore, there are 9 transitions in total, corresponding to all possible arrows of Figure [T| We 
list them below, stressing that synchronization between clients and servers has a rate computed 
using the minimum, in the PEPA style [3lj. With Ix we denote a vector of length n which is 
equal to 1 for component X and zero elsewhere. 

• request, iwequest — X^rq ^ L^ui; ^rq ^ '-'pj > "^request — "^^ ' rnml^Cyj-g, O^qj).! 

• reply: Kreply ^ \^w ^ ^tjJrp ^ Sl\, Vreply = T^^T^ykw^wj '^rp^rp)', 

• timeout (client): Rumeoutl = {Cw — ^ Crc}, fUmeoutl = hoCu,', 

• recover: itrecover -- X^rc ' ^rqi^ "^recover ~~ f^rec^rc] 

• think: Rthink = {Ct ^ Crq, Srp -^ Si}, r think = ktCt, 

• logging: Rlogging = {Si -^ Srq}, Rlogging = hSl] 

• process. 'V process — X'-'p ' '-'rpj j ^process — "^p^pj 

• timeout (server processing): Rtimeout2 = {Sp -^ Si}, rtimeout2 = kstoSp, 

• timeout (server replying): Rtimeouts = {Srp^ Si}, mmeouts = kstoSrp] 



The system-level models we have defined depend on the total population A^ and on the ra- 
tion between server and clients, which is specified by the initial conditions. Increasing the total 
population A^ (keeping fixed the client-server ratio), we obtain a sequence of models, and we are 
interested in their limit behaviour, for A^ going to infinity. 

In order to compare the models of such a sequence, we will normalize them to the same scale, 
dividing each variable by A^ and thus introducing the normalized variables X^-^^ = jy . In the case 
of a constant population, normalised variables are usually referred to as the occupancy measure, 
as they represent the fraction of agents in each state. Update vectors are scaled correspondingly, 
i.e. dividing them by A^. Furthermore, we will also require a proper scaling (in the limit) of the 
rate functions of the normalized models. More precisely, let X^^' = (X^^',?"*-^', Xq' ') be the 
A^-th non-normalized model and X'-"' = (X' ^T' %Xq ) the corresponding normalized model. 
We require that: 

• initial conditions scale appropriately: Xq = —^ — ; 

• for each transition {\r^,rr (X)) of the non-normalized model, we let fr (X) be the rate 
function expressed in the normalised variables (i.e. after a change of variables). The corre- 
sponding transition in the normalized model is {RT^rr (^)), with update vector equal to 
jyV,-. We assume that there exists a bounded and Lipschitz continuous function /r(X) : 
E — )• M" on normalized variables (where E contains all domains of all X^ >), independent of 

A^, such that ^^ ^^' — )• /t(x) uniformly on E. 

In accordance with the previous subsection, we will denote the state of the CTMC of the A^-th 
non-normalized (resp. normalized) model at time t as 'K^'^'{t) (resp. 'K^^'{t)). 

Example. Consider again the running example. If we want to scale the model with respect to the 
scaling parameter A'^, we can increase the initial population of clients and servers by a factor k 
(hence keeping the client-server ratio constant), similarly to [35j. The condition on rates, in this 
case, automatically holds due to their (piecewise) linear nature. For non-linear rate functions, the 
convergence of rates can usually be enforced by properly scaling parameters with respect to the 
total population A^. 

3.2 Deterministic limit theorem 

In order to present the "classic" deterministic limit theorem, we need to introduce a few more 
concepts needed to construct the limit ODE. Consider a sequence of normalized models X^ ' and 
let v,- be the (non- normalised) update vectors. The drift F^^' (X) of X is defined as 

F^''\y^) = Y.l^-rr^r''\y^) (2) 

Furthermore, let fr'-E^ M", r G T" be the limit rate functions of transitions of X^^' . We define 
the limit drift of the model X^^' as 

F(X) = J]v.A(X) (3) 

ret 

It is easily seen that ^(^^(x) — )• -F(x) uniformly. 

The hmit ODE is ^ = F(x), with x(0) = xq G 5. Given that F is Lipschitz in E (as ah fr 
are), the ODE has a unique solution x(t) in E starting from xq. Then, the following theorem can 
be proved [71 [H] : 



Theorem 3.1 (Deterministic approximation [7, 8j). Let the sequence ^^^'{t) of Markov processes 
and x(t) be defined as before, and assume that there is some point xq G 5" such that X*-^'*(0) — )■ xq 
in probability. Then, for any finite time horizon T < oo, it holds that: 

f\ sup ||x(^)(i)-x(t)|| >ei ^0. 

\^0<t<T J 

Notice that the theorem can be speciahsed to subsets E' (^ E, in which case it can also provide 
an estimate of exit times from set E' , see [8]. Furthermore, if the initial conditions converge almost 
surely, then it also holds that supo<i<T 11^ (^) ~ x(f)|| — ;• almost surely [36] , 

3.3 Fast simulation 

We now turn our attention back to a single individual in the population. Even if the system-level 
dynamics, in the limit of a large population, becomes deterministic, the dynamics of a single agent 
remains a stochastic process. However, the fluid limit theorem implies that the dynamics of a 
single agent, in the limit, becomes essentially dependent on the other agents only through the 
global system state. This asymptotic decoupling allows us to find a simpler Markov Chain for the 
evolution of the single agent. This result is often known in the literature [9j under the name of 
fast simulation 



To explain this point formally, let us focus on a single individual Y^ , which is a Markov 
process on the state space S = {l,...,n}, conditional on the global state of the population 
X^ '(t). Let Q' •'(x) be the infinitesimal generator matrix of Y^ , described as a function of the 
normalized state of the population X'^ = x, i.e. 

p{y,(^)(t + dt) = j I yi^)(t) = i, xW(t) = x} = qif{^)dt. 

We stress that this is the exact Markov Chain for Y"^ , conditional on X*^ ■'(t), and that this 
process is not independent of 'K^^'{t). In fact, without conditioning on X*-^', Y"^ (t) is not a 
Markov process. This means that in order to capture its evolution in a Markovian setting, one has 
to consider the Markov chain (Y^ (t), X' •'(t)). 

Example. Consider the running example, and suppose we want to construct the CTMC for a single 
client. For this purpose, we have to extract from the specification of global transitions a set of local 
transitions for the client. The state space of a client will consist of four states, Sc = {rq,w,t,rc}. 

Then, we need to define its rate matrix Q^ ' . In order to do this, we need to take into account 
all global transitions involving a client, and then extract the rate at which a specific client can 
perform such a transition. As a first example, consider the think transition, changing the state of 
a client from t to rq. Its global rate is rthink = ktCt- As we have Ct clients in state t, the rate at 
which a specific one will perform a think transition is -^ = kt- Hence, we just need to divide the 
global rate of observing a think transition by the total number of clients in state t. Notice that, 
as we are assuming that one specific client is in state i, then Ct > 1, hence we are not dividing by 
zero. 

Consider now a reply transition. In this case, the transition involves a server and a client in 
state w. The global rate is rrepiy = min{kwCw, krpSrp), and C^ > 1 (in the non- normalized model 

C 

with total population N). Dividing this rate by Cw, we obtain min(fc^, fcrp^), which is defined 
for Cw > 0. If we switch to normalised variables, we obtain a similar expression: min(fc^, fc^p-^), 
which is independent of A^. However, in taking N to the limit we must be careful: even if in the 



non-normalized model C^ (and hence c^) are always non-zero (if a specific agent is in state w), 
this may not be true in the limit: if only one client is in state w, then the limit fraction of clients 
in state w is zero (just take the limit of -^). Hence, we need to take care of boundary conditions, 
guaranteeing that the single-agent rate is defined also in these circumstances. In this case, we can 
assume that the rate is zero if Srp is zero (whatever the value of c„,), and that the rate is k^ if c^ 
is zero but Srp > 0. 

In order to treat the previous set of cases in a homogeneous way, we make the following 
assumption about rates: 

Definition 3.1. Let r G 7" be a transition such that its update rule set contains the rule i -^ j, 
with multiplicity nir^i-^j- The rate r,- is single- agent-i compatible if there exists a Lipschitz 
continuous function /^(x) on normalized variables such that the limit rate on normalized variables 
/t-(x) can be factorised as /t(x) = Xj/^(x). A transition r is single- agent compatible if and only if 
it is single- agent-i compatible for any i appearing in the left-hand side of an update rule. 

Hence, the limit rate of observing a transition from i to j for a specific agent in state i is 
mr,i^jfl{x), where the factor rur^i-^j comes from the fact that it is one out of rur^i^j agents 
changing state from i to j due to rjj 

Then, assuming all transitions r are single-agent compatible, we can define the rate q\ ■ as 

tGT I {i^jKRr * ret I {i^jjCRr * 

It is then easy to check that 






reT I {i^j}CB.r * ret I {i^jjCRr 

In the following, we fix an integer /c > and let Z^^^ = (yI^\ . . . ,Y^'^^) be the CTMC 
tracking the state of k selected agents among the population, with state space S = S^. Notice 
that k is fixed and independent of A^, so that we will track k individuals embedded in a population 
that can be very large. 

Let x(t) be the solution of the fluid ODE, and assume to be under the hypothesis of Theorem 
Consider now zj^ (t) and Zk(t), the time-inhomogeneous CTMCs on S defined by the following 



3.1 



infinitesimal generators (for any h = 1, . . . ,k): 

P{zf)(t + dt) = izu...,j,...,zk)\ zf )(t) = (zi, . . . , i, . . . , Zk)} = gjf (x(i))dt, 

¥{zk{t-^dt) = {zi,...,j,...,Zk) I Zk{t) = {zi,...,i,...,Zk)} = qij{x{t))dt, 

Notice that, while Zj^ describes exactly the evolution of k agents, z^. and z^ do not. In fact, 
they are CTMCs in which the k agents evolve independently, each one with the same infinitesimal 
generator, depending on the global state of the system via the fluid limit. 
However, the following theorem can be proved [9j: 

Theorem 3.2 (Fast simulation theorem). For any T < oo, ^{Z^ (t) 7^ zj, (t), for some t < 
T} -^ 0, and F{Zk{t) / Zk{t), for some t<T} ^0, as N ^ 00. 

^The factor m stems from the following simple probabilistic argument: if we choose at random m agents out of 
Xi , then the probability to select a specific agent is ^ . 



This theorem states that, in the hmit of an infinite population, each fixed set of k agents wih 
behave independently, sensing only the mean state of the global system, described by the fiuid limit 
x(t). Furthermore, those k agents will evolve independently, as if there was no synchronisation 
between them. This asymptotic decoupling of the system, holding for any set of k agents, is also 
known in the literature under the name of propagation of chaos [6J. In particular, this holds if 
we define the rate of the limit CTMC either by the single-agent rates for population A^ (z^ ) 
or by the limit rates {zk). Note that, when the CTMC has density dependent rates [36j, then 
zj^ [t) = Zk{t), as their infinitesimal generators will be the same. 

We stress once again that the process Zj^ (t) is not a Markov process. It becomes a Markov 
process when considered together with X*^ •'(t). This can be properly understood by observing 
that it is the projection of the Markov process [Y^ (t), . . . , Y^ (t)) on the first k coordinates, 
and recalling that a projection of a Markov process need not be Markov (intuitively, we can throw 
away some relevant information about the state of the process). However, being the projection 
of a Markov process, the probability of Z^. (t) at each time t is perfectly defined. Nevertheless, 
its non-Markovian nature has consequences for what concerns reachability probabilities and the 
satisfiability of CSL formulae. 

Example. Consider again the client-server example, and focus on a single client. As said before, 
its state space is Sc = {rq, w, t, re}, and the non-null rates of the infinitesimal generator Q for the 
process zi are: 

• qrq,w{t) = kr ram{l, Srq{t) / Crq{t)} (with appropriate boundary conditions); 

• Qw,t{t) =mm{kw,krpSrp{t)/cu,{t)}; 

• Qw,rc{j^) — "^to: 

• Qt,rq{t) = h; 

• Qrc^rqxi^) -- farc- 
in Figure [2j we show a comparison of the transient probabilities for the approximating chain 

for a single client and the true transient probabilities, estimated by Monte Carlo sampling of the 
CTMC, for different population levels A^. As we can see, the approximation is quite precise already 
for A^ = 15. 

Remark 3.1. Single-agent consistency is not a very restrictive condition. However, there are cases 
in which it is not satisfied. One example is passive rates in PEPA |3l]. In this case, in fact, the 
rate of the synchronization of P = (a,T).Pl and Q = {a,r).Ql is rXQl{Xp > 0}. In particular, 
the rate is independent of the exact number of P agents. If we look at a single P-agent rate, it 
equals rj^l{Xp > 0}. Normalising variables, we get the rate r-^l{xp > 0}, which approaches 
infinity as xp goes to zero (for xq fixed). Hence, it cannot be extended to a Lipschitz continuous 
function. However, in the case xp = and xq > 0, if we look at a single agent, then the speed at 
which P changes state is in fact infinite. We can see this by letting Xp = 1 and Xq = Nq, so that 
the rate of the transition from the point of view of P is Nq -^ oo. Thus, in the limit, the state P 
becomes vanishing. 

Remark 3.2. The hypothesis of constant population, i.e. the absence of birth and death, can 
be relaxed. The fiuid approximation continues to work also in the presence of birth and death 
events, and so does the fast simulation theorem. In our framework, birth and death can be easily 
introduced by allowing rules of the form — >• i (for birth) and i — ?■ (for death). In terms of a 
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Figure 2: Comparison of tlie transient probability for all four states of the fluid model of the client 
server system, computed solving the Kolmogorov forward equations, and the transient probability 
of CTMC models for A^ = 15 and A^ = 150 (2:1 client server ratio). Parameters are kr = 1, 



/C), 



100, kto = 0.01, kt = 1, krc = 100 ki = 10, kp = 0.1, k. 



sto 



0.005, initial conditions of the 



full system are Crq = n, Srq = m, while the single client CTMC starts in state rq. 
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single agent, death can be dealt with by adding a single absorbing state to its local state space 
S. Birth, instead, means that we can choose the time instant at which an agent enters the system 
(provided that there is a non-null rate for birth transitions at the chosen time). 

Another solution would be to assume an infinite pool of agents, among which only finitely 
many can be alive, and the others are an infinite supply of "fresh souls". Even if this is plausible 
from the point of view of a global model, it creates problems in terms of a single agent perspective 
(what is the rate of birth of a soul?). A solution can be to assume a large but finite pool of agents. 
But in this case birth becomes a passive action (and it introduces discontinuities in the model, 
even if in many cases one can guarantee to remain far away the discontinuous boundary), hence 
we face the same issues discussed in Remark 13.11 



3.4 Continuous Stochastic Logic 

In this section we consider labelled stochastic processes. A labelled stochastic process is a random 
process Z{t), with state space S and a labelling function L : 5 — t- 2^, associating with each state 
s £ S a subset of atomic propositions L{s) C V = {ai, . . . , a^ . . .} true in that state: each atomic 
proposition Oj S "P is true in s if and only if aj G L{s). We require that all subsets of paths 
considered are measurable. This condition will be satisfied by all subsets considered in the paper, 
as Z{t) will always be either a CTMC, defined by an infinitesimal generator matrix Q{t) (possibly 
depending on time), or the projection of a CTMC. From now on, we always assume we are working 
with labelled stochastic processes. 

A path of Z{t) is a sequence a = sq — ^ si -^ . . ., such that the probability of going from Si 
to Sj+i at time tfj[i] = '}2)=o'tj^ i^ greater than zero. For CTMCs, this condition is equivalent to 
Qsi,Si+i{ta[i]) > 0. Denote with a@t the state of a at time t, with a[i] the i-th state of a, and with 
tij[i] the time of the i-th jump in a. 

A time-bounded CSL formula 93 is defined by the following syntax: 



"P 



a I (^1 A v?2 I -V' I Pmp(X[^1'^2V) I n<p(¥'iU[^i'^2l(^2). 



The satisfiability relation of <p with respect to a labelled stochastic process Z{t) is given by the 
following rules: 



s,to 
s,to 
s,to 
s,to 
s,to 



a if and only if a G L{s); 

-193 if and only if s, to ^ ip; 

If I A ip2 if and only if s, to |= ^1 and s, to |= 9^2; 

Pmp(X[^1'^2](^) if and only if¥{a \ a, to |= Xl^i'^^V} ^P- 

i^Mp(viU['^i'^2]y,2) if and only if Pjcj | a, to |= ipilJ^^^'^^^ip2} ^p. 



a, to 1= Xl^i'^^V if and only if t<^[l] € [Ti, Ta] and f7[l], to + t^[l] 



if. 



• cr,to ^ v3iU['^i'^2](^2 if and only if 3t G [to + ri,to + T2] s.t. cr@t,t |= (^2 and Vto < t < t, 
a@t,t 1= (fi. 

Notice that we are considering a fragment of CSL without the steady state operator and allow- 
ing only time-bounded properties. This last restriction is connected with the nature of convergence 



theorems 3.1 and 3.2, which hold only on finite time horizons. However, see Remark 6.3 for possible 
relaxations of this restriction. 
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Model checking of a next CSL formula Pt^ipi^^'^'^''^''^) is usually performed by computing the 
next-state probability via an integral, and then comparing the so-obtained value with the threshold 
p. Model checking of an until CSL formula Ptxip{y^i'U^'^^''^^'ip2) in a time-homogeneous CTMC Z{t), 
instead, can be reduced to the computation of two reachability problems, which themselves can 
be solved by transient analysis [i3|. In particular, consider the sets of states U = l^^pi} and 
G = lip2l and compute the probability of going from state si [/ to a state S2 ^ U in Ti 
time units, in the CTMC in which all [/-states are made absorbing, vr^^^2(Ti). Furthermore, 
consider the modified CTMC in which all U and G states are made absorbing, and denote by 
TTgj, 53(^2 — Ti) the probability of going from a state S2 ^ U to a state S3 G G in T2 — Ti units of 
time in such a CTMC. Then the probability of the until formula in state s can be computed as 
Psiv) = X^sasG S2^c/'''"si, 52(^1)^52, S3 (^2 — ^i)- The probabilities ir^ and tt^ can be computed using 
standard methods for transient analysis (e.g. by uniformisation [37] or by solving the Kolmogorov 
equations [38]). Then, to determine the truth value of the formula if in state s, one has just to solve 
the inequality Ps{f) >^ P- The truth value of a generic CSL formula can therefore be computed 
recursively on the structure of the formula. 

4 Next State Probability 

In this section, we start the presentation of the algorithmic procedures that underlie the CSL 
model checking algorithm. We will focus here on next state probabilities for single agents (or a 
fixed set of agents), in a population of growing size. In particular, we will show how to compute 
the probability that the next state in which the agents jumps belongs to a given set of goal states 
G C 5, constraining the jump to happen between time [to + ^1,^0 + T2], where to is the current 
time. This is clearly at the basis of the computation of the probability of next path formulae. 
More specifically, we provide algorithms for ICTMC (hence for the limit model Zk{t)), focussing 
on two versions of the next state probability: the case in which the set G is constant, and the case 
in which the set G depends on time (i.e. a state may belong to G or U depending on time t). 

Definition 4.1. Let Z{t) be a CTMC with state space S and infinitesimal generator matrix Q{t). 

1. Let G C S. The constant-set next state probability Pnext{Z,to,Ti,T2, G)[s] is the probability 
of the set of trajectories of Z jumping into a state in G, starting at time to in state s, within 
time [to + Ti,tQ + T2]. Pnext{Z, to, Ti, T2, G) is the next state probability vector on S. 

2. Let G : [to,ti] x S — ;■ {0,1} be a time-dependent set, identified with its indicator function 
(i.e. G{t is the goal set at time t). 

The tim,e-varying-set next state probability PnextiZ,to,Ti,T2,G{t))[s] is the probability of 
the set of trajectories of Z jumping into a state in G{t) at time t S [to -|- Ti , to -|- T2], starting 
at time to in state s. 

The interest in the time-varying sets is intimately connected with CSL model checking. In 
fact, the truth value of a CSL formula in a state s for the (non-Markov) stochastic process Zf, (t) 
depends on the initial time at which we start evaluating the formula. This is because Zj^ (t) 
depends on time via the state of X' •'(t). Furthermore, time- dependence of truth values of 
CSL formulae manifests also for the limit process Zk{t), which is a time-inhomogeneous Markov 
Process. Therefore, in the computation of the probability of the next operator, which can be 
tackled with the methods of this section, we need to consider time-varying sets. As a matter of 
fact, we will see that dealing with time- varying sets (like those obtained by solving the inequality 
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PnextiZ, to, Ti, T2, G{t)) ixi p, for CXI G {<, <, >>}) requires us to impose some additional regularity 
conditions on the rate functions of Z and on the time-dependence of goal sets. 

In the following sections, we will first deal with the computation of next state probabilities for 
a generic ICTMC Z{t) and then study the relationship between those probabilities for Zj^ (t) and 

Zkit). 

4.1 Computing next-state probability 

Consider a generic ICTMC Z{t), and focus for the moment on a constant set G C 5. For any fixed 
to, the probability Pnext{Z, tQ,Ti,T2, G)[s] that Z(t)'s next jump happens at time t £ [to + 7i, to + 
T2] and ends in a state of G, given that Z{t) is in state s at time to, is given by the following 
integral [Ml [25] 



PnextiZ,to,Ti,T2,G)[s]= (z.,G(t) • e-^(*«'*)Wdt, (4) 

Jto+Ti 

where A(to,t)[s] = / —qs,s(,T)dT, is the cumulative exit rate of state s from time to to time t, and 

Jto 
Qs,G{t) = Yls'eG^s,s'{t) is the rate of jumping from s to a state s' £ G at time t. 

Equation [4| holds for the following reason. Let At be the event that we jump into a G state 

in a time r G [to,t]. Then At^ C At^ for h < ta, and F{At} = / qs,G{t) • e"^(*0'*)Wdt. We are 

Jto 

interested in the probability of the event A = ^to+T2 \ ^to+Ti, which has probability 
F{A} = nAto+T,} - IP{A,+tJ = r^ ' qs,G{t) ■ e-^(*°'*)Wdt. 

Jto+Ti 

In order to compute PnextiZ,tQ,Ti,T2,G)[s] for a given to, we can numerically compute the 
integral, or transform it into a differential equation, and integrate the so-obtained ODE with 
standard numerical methods. This simplifies the treatment of the nested integral A(to, t)[s] involved 
in the computation of Pnext- More specifically, we can introduce two variables, P and L, initialise 
P(to + Ti) = and L(to + Ti) = A(to,to + Ti), and then integrate the following two ODEs from 
time to + Ti to time to + T2 : 

^Pit) = qsMt) ■ e-^(*) 

d (^) 

-L{t) = -qsAt) 

However, for CSL model checking purposes, we need to compute 
Pnext{Z,to,Ti,T2,G)[s] as a function of to : -Ps(to) = Pnext{Z,to,Ti,T2,G)[s]. One way of doing this 
is to compute the integral Q for any to. A better approach is to use the differential formulation 
of the problem, and define a set of ODEs with the initial time to as independent variable. First, 
observe that the derivative of -Ps(to) with respect to to is 



dt, 



-^s(to) = qs,G{t0 + T2) ■ e-^^'o,to+T,) _ ^^^^(^^ ^ ^^) . g-A(to,to+Ti) 







rto+T2 f) . , , 

+ / 7^g.,G(t)-e-^(*-*)dt 

JtQ+Ti OtQ 
qs,G{to + T2) ■ e-A(*0'*o+^^) - g,,G(to + Ti) • e-^(*0'*o+^i) 
- qs,s{to)Ps{to) 
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function next-state-probability(Z, G, Ti, T2, to, ^i) 
for all s E 5 do 

Compute Ps(to) by solving ODEJs] 
Compute Ps{t) for t £ [to,ti] by solving ODEJg] 
end for 
return P(t), t £ [to,ti] 
end function 

Figure 3: Algorithm for the computation of next-state probability P{t), for any state s and 
t G [fo,ti]- Other input parameters are as in the text. 

Consequently, we can compute the next-state probability as a function of to by solving the following 
set of ODEs: 

|n(t) = qs,G{t + T2) • e-^2W _ g^_^(i + Ti) • e-^iW - g.,.(t)P,(t) 

|Li(t) = -qs,s{t) + QsAt + T^i) (6) 

^L2(t) = -g,,,(t) + g,,s(t + r2) 

where Li(t) = A(t,t + Ti) and L2(t) = A(t,t + r2). 

Initial conditions are Ps(to) = PnextiZ,to,Ti,T2,G)[s], -Li(to) = A(to,to + Ti), and L2(to) = 

A(to,to + T2), and are computed solving the equations [5j The algorithm is sketched in Figure [SJ 

We turn now to discuss the case of a time- varying next-state set G(t). In this case, the only 
difference with respect to the constant-set case is that the function q.^ft) is piecewise continuous, 
rather than continuous. In fact, each time a state s' gains or loses membership of G{t), the range 
of the sum defining q.^G{t) changes, and a discontinuity can be introduced. However, as long as 
these discontinuities constitute a set of measure zero (for instance, they are finite in number), this 
is not a problem: the integral Q is defined and absolutely continuous, and so is the solution of 
the set of ODEs ([6| (because the functions involved are discontinuous with respect to time). It 
follows that the method for computing the next-state probability for constant sets works also for 
time- varying sets. 

Now, if we want to use this algorithm in a model checking routine, we need to be able also to 
solve the equation Ps{t) = p, for p £ [0, 1] and each s G 5. In particular, for obvious computability 
reasons, we want the number of solutions to this equation to be finite. This is unfortunately not 
true in general, as even a smooth function can be equal to zero on an uncountable and nowhere 
dense set of Lebesgue measure (for instance, on the Cantor set [40j ) . 

Consequently, we have to introduce some restrictions on the class of functions that we can use. 
In particular, we will require that the rate functions of z^ and of Zj^ are piecewise real analytic 
functions. 

4.1.1 Piecewise Real analytic functions 

A function / : / — t- M, / an open subset of M, is said to be analytic [H] in / if and only if for 
each point to of / there is an open neighbourhood of / in which / coincides with its Taylor series 
expansion around to- Hence, / is locally a power series. For a piecewise analytic function, we 
intend a function from / — t- M, / interval, such that there exists Ii, . . . ,Ik disjoint open intervals, 
with I = [jj Ij, such that / is analytic in each /-,-. A similar definition holds for functions from M" 
to R, considering their multi-dimensional Taylor expansion. 
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Analytic functions are a class of functions closed by addition, product, composition, division 
(for non-zero analytic functions), differentiation and integration. Piecewise analytic functions also 
satisfy these closure properties, by considering the intersections of their analytic sub-domains. 
Many functions are analytic: polynomials, the exponential, logarithm, sine, cosine. Using the 
previous closure properties, one can show that most of the functions we work with in practice are 
analytic. 

Analytic functions have two additional properties that make them particularly suitable in this 
context: 

1. The zeros of an analytic function / in /, different from the constant function zero, are 
isolated. In particular, if / is bounded, then the number of zeros is finite. This is true also 
for the derivatives of any order of the function /. 



di 

consequence of the Cauchy-Kowalevski theorem [42J 



2. If / is analytic in a set E, then the solution x of -^ = /(x) in E is also analytic (this is a 



This second property, in particular, guarantees that if the rate functions of Zk and Z^ are 
piecewise analytic, then all the probability functions computed solving the differential equations, 
like those introduced in the previous section, are also piecewise analytic. 

In the following, we will need the following straightforward property of piecewise analytic 
functions: 

Proposition 4.1. Let f : I ^M be a piecewise analytic function, with / C M a compact interval. 
Let Sj = {x G M I ^^(/~^({x})) = 0} he the set of all values x such that f is not locally constantly 
equal to x, where fi£ is the Lebesgue measure. Furthermore, let Z^ = f~^{{x}) be the set of 
solutions of f{t) = X and let DZf = {x € M | Vt G Zx, fit) ^ 0}. Then 

1. \/x £ Ef, Zx is finite. 

2. fii{Ef n DZf) = 1 

4.2 Convergence of next-state probability 

We consider now the problem of relating the next-state probabilities for the limit single agent 
process Zk{t) and the sequence of single agent processes Zj^ (t) in a population of size A^. In 
particular, we want to show that the probability Pg (t) = Pnext{Z\. \tQ^Ti,T2,G)[s\ converges 
to Ps{t) = Pnext{zk-,to-,Ti,T2,G)[s\ Uniformly for t G [tojii]; as A^ goes to infinity. We will prove 
this result in a general setting. More specifically, we will consider time- varying sets that can depend 
on A^, and that converge to a limit time-varying set in a suitable sense. This is needed because 
the time-varying sets we need to consider are obtained by solving (for each s G 5) equations of 
the form Pg (t) — p = or Ps{t) — p = 0, which are generally different, but intuitively converge 
(as Ps (t) converges to Ps{t)). 

4.2.1 Robust time-varying sets 

We first introduce a notion of robustness for time-varying sets, which will be enforced on limit 
sets: 

Definition 4.2. A time-dependent subset V{t) of 5, t G /, is robust if and only if there is a 
piecewise analytic function /ly : 5 x / — )• M and an operator txiG {<,<,>,>}, such that for each 
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s £ S, the indicator function V^ : / — )■ {0,1} of s is given by l{/iy(s,t) M 0}, and it further 
satisfies: 

1. the number of discontinuity points Disc{V) = {{s,t} \ Vs{t~ ^ Vs{i^)} is finite; 

2. if hv is analytic in {s,t) and hy{s,t) = 0, then ^hv{s,t) / (zeros of hy are simple); 

3. if hy is not analytic in {s,t), then hv{s,t^) / and hv{s,t~^) / 0|j 

In the following, we will usually indicate with V{t) both a time dependent set V and its 
indicator function (with values in {0, 1}™, m = \S\), and with hy the piecewise analytic function 
defining it. 

As we will see later on, the notion of robustness is closely related to the comput ability of the 
reachability probability for time-varying sets and with the decidability of the model checking 
algorithm for ICTMCs, both discussed in Section [6| 

Furthermore, we need the following notion of convergence for time-varying sets: 

Definition 4.3. A sequence of time-varying sets V^ '(t), t £ I, converges robustly to a robust 
time- varying set V{t), t S I, if and only if, for each s £ S and each open neighbourhood U of 
Disc{Vs) (i.e. the set of discontinuity points of Vg), V} (t) — )• Vs{t) uniformly in I \U. 



Connecting the notions of robust set and robust convergence, we have the following: 

Proposition 4.2. Let V^ '{t) be a sequence of time varying sets converging robustly to a robust 
V 



set V{t), t G /. Let L»[f ^ = {t \ V^^\t) / V{t)}. Then fie{D'y'^) -^ 0, where /X£ is the Lebesgue 



measure on 



4.2.2 Convergence results 

We are now ready to state the convergence result for next-state probabilities. We will assume that 
the limit time- varying set is robust. The following lemma will be one of the key ingredients to 



prove the inductive step in the convergence for truth of CSL formulae, see Section 6.2 

Lemma 4.1. Let X^ ' be a sequence of C TMC models, as defined in Section 3.1. and let Z\, 

and Zk be defined from X^^' as in Section 3.3, with piecewise real analytic rates, in a compact 

interval [0, T'], for T' >ti+Tb. 

Let G{t), t G [to,ti + Ti,] be a robust time-varying set, and let G^^'{t) be a sequence of time-varying 

sets converging robustly to G. 

Furthermore, let P{t) = Pnextizk,t,Ta,Tf,,G) and 

pW{t) = PneAzi''\t,Ta,n,G), t G [t^M]- 

Finally, fix p & [0,1], OOG {<,<,>,>}, and let Vp{t) = l{P{t) Mp}, vi^\t) = l{pW(t) Mp}. 
Then 

1. P^^\t) -^ P{t), uniformly in t G [to,ti]. 

2. For almost every p G [0, 1], Vp is robust and the sequence Vp converges robustly to V^, 



p- 



^This condition states that, if hv is continuous but not analytic in (s,i), then it cannot be equal to zero in 
those points, implying that first order derivatives exist and are non-null in all continuity points in which hv crosses 
zero. Moreover, if hv{s,t~) ^ hviSjf^), hv can cross zero in (s,t) only if the jump contains zero, meaning that 
mm{hv{s,t~),hv{s,t'^)} < < ma,x{hv{s,t~),hv{s,t~^)} . 
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5 Reachability 

In this section, we will focus on the computation of reachability probabilities of a single agent 
(or a fixed set of agents), in a population of increasing size. Essentially, we want to compute the 
probability of the set of traces reaching some goal state G C 5 within T units of time, starting at 
time to and avoiding unsafe states in [/ C 5. The key point is that the reachability probability 



of the limit CTMC Zk{t) obtained by Theorem 3.2 approximates the reachability probability of a 



single agent in a large population of size A^, i.e. the reachability probability for Zj^ (t). 

Similarly to Section [4j we will consider two versions of the reachability problem: one for 
constant goal and unsafe sets, and one in which G and U depend on time. We will state these 
problems for a generic ICTMC Z{t) on state space S: 

Definition 5.1. Let Z{t) be an ICTMC with state space S and infinitesimal generator matrix 
Q{t). 

1. Let U,G CS. The constant-set reachability Preachi^, to,T, G , U)[s] is the probability of the 
set of trajectories of Z reaching a state in G without passing through a state in U, within T 
time units, starting at time to in state s. Preachi^, to, T, G, U) is the reachability probability 
vector on S. 

2. Let U, G : [to, ti] x 5 — )• {0, 1} be time-dependent sets, identified with their indicator function 
(i.e. G{t), U{t) are the goal and the unsafe sets at time t). The time-varying-set reachability 
PreachiZ, to, T, G{t), U{t))[s] is the probability of the set of trajectories of Z reaching a state 
in G{t) at time t € [to, to + T] without passing through a state in U{t'), for t' G [to, t], starting 
at time to in state s. 

In the following sections, we will first deal with the specific reachability problem for a generic 
ICTMC Z{t), presenting an effective way of computing such probability, and then studying the 
relationship between the reachability probabilities of Zj^ (t) and Zk{t)- 

5.1 Constant-set reachability 

We consider constant-set reachability, according to Definition |5.1[ For the rest of this section let 
Z{t) be an ICTMC on 5, with rate matrix Q{t) and initial state Z{Qi) = Zo £ S. We will solve the 
reachability problem in a standard way, by reducing it to the computation of transient probabilities 
in a modified ICTMC |13j . The solution is similar to the one proposed in [26j. 

Let n(ti,t2) be the probability matrix of Z(t), in which entry 715^^52(^1)^2) gives the probability 
of being in state S2 at time t2, given that we were in state si at time ti. The Kolmogorov forward and 
backward equations describe the time evolution of n(ti, t2) as a function of t2 and ti, respectively. 
More precisely, the forward equation is — q^' = n(ti, t2)(5(i2), while the backward equation is 

^%*^ = -Q(ti)n(ti,t2). 

The constant-set reachability problem, for a given initial time to, can be solved by integra- 
tion of the forward Kolmogorov equation (with initial value given by the identity matrix) in the 
modified ICTMC Z'(t), with infinitesimal generator matrix Q'{t), in which all unsafe states and 
goal states are made absorbing [13j (i.e. q's^^ S2^^) ~ ^' ^^^ each si £ G U U). In particular, 
Preach{Z,to,T,G,U) = n'(to,to + T)eG, where e^ is an n x 1 vector equal to 1 if s G G and 
otherwise, and 11' is the probability matrix of the modified ICTMC Z'£j We emphasise that, in 

^Clearly, alternative ways of computing the transient probability, like uniformization for ICTMC [43], could also 
be used. However, we stick to the ODE formulation in order to deal with dependency on the initial time to- 
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function reachability-constant-set(Z, T, G, U, to, h) 

Construct the CTMC in which G and U states are absorbing, with rate matrix Q'{t). 

Compute n'(to,io + T) by solving the forward Kolmogorov ODE for the modified CTMC. 

Compute n'(t,f + T) for t € [to,ii] by solving ODE for the modified CTMC with initial 
conditions n'(to, to + T). 

return P{t) = n'(io, to + T)eG, t G [io, ti]. 
end function 

Figure 4: Algorithm for the computation of reachability probability P{t) for t G [toj ^i] and constant 
goal and unsafe sets G and U. Other input parameters are as in the text. 

order for the initial value problem defined by the Kolmogorov forward equation to be well posed, 
the infinitesimal generator matrix Q{t) has to be sufficiently regular (e.g. bounded and integrable). 
As already remarked, in contrast with time- homogeneous CTMC, the reachability probability 
for ICTMC can depend on the initial time to at which we start the process. Consider now the 
problem of computing P{t) = PreachiZ, t, T, G, U) as a function of t G [to, ti]. To this end, we can 
solve the forward equation for to aiid then use the chain rule to define a differential equation for 
n(t, t + T), solving it using n(to, to + T) as the initial condition, i.e. 

dU{t, t + T) _ an(t, t + T) dli{t, t + T) d{t + T) 

Jt ~ dt ^ d{t + T) Jt (7) 

= -Q(t)n(t, t + T)+ U{t, t + T)Q{t + T). 

Using a numerical solver for the ODE, this gives an effective algorithm (Figure |4| to compute 
the probability of interest (for any fixed error bound). Furthermore, if we can guarantee that the 
number of zeros of the equation P{t) — p is finite, then we also have an effective procedure to 
compute the truth value of P{t) to p, for cxiG {<, <, >, >} (provided we can find those zeros, as 
will be discussed in Section |6l). 



Consider now the sequence of processes Z^ defined in Section 



3.3, We are interested in the 



asymptotic behaviour of Preach{Z), ,t,T, G, U). The following result is a consequence of Theorem 
3^1 

Proposition 5.1. Let X^"' he a sequence o f CT MC models, as defined in Section 



3.1 



and let Z\, 



and Zk be defined from X^^' as in Section 3.3. Assume that the infinitesimal generator matrix 



Q{t) of Zk is bounded and integrable in every compact interval [0,T]. Then 

Preach{z[ ,t,T, G, U) -> Preach{zk,t, T, G, U), Uniformly in [to, ti], as N -j- oo 

i.e. SUPtg[tj,_t^] \\PreachiZk ,t,T,G,U) - PreachiZk,t,T,G,U)\\ -^ 0. 

The previous proposition shows that the reachability probability for Zj^ converges to the 
reachability probability for z^, hence for large N we can approximate the former with the latter. 
It is interesting to observe how the reachability probability for Z^ (t) depends on the initial 



time. As previously remarked, Zj^ (t) is not a Markov- process, but {Zj^ ,X'-^')(t) is. Fur 

^k (*) ^y projecting on the first component of (Z^ 



thermore, we can obtain Zi (t) by projecting on the first component of (Zi ,X'^-')(t). The 



reachability probability for Z^ (t) can be obtained from that of (Z^ ,'K^ '){t) in the follow- 
ing way: compute the reachability probability Pu^g{s,x,T) for each state (s,x) of (Z^ ,X(^)) 
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with time horizon T. As {Z^ ,'K^^') is a time-homogeneous CTMC, this probabiUty is inde- 
pendent of the initial time. Fix a state s G 5 of Z^ , and consider the probabihty -Ps,x(i|s) = 
P{(Z^ ,X^ ')(t) = (s, x) I Zj^ (t) = s} of being in (s,x) at time t, conditional on being in s, 
i.e. Ps,x(t|s) = P{(zi^\x(^))(t) = {s,x)}/J2^F{{zl^\±^^^){t) = {s,x)} (when the denomina- 
tor is non-zero). Then, this is the initial distribution of (Z^ ,'K^^') that we have to take into 

account when computing the reachability probability Pj.each{Z\. ,t,T, G, U)[s\, starting at time t. 
It follows that 

Preach{Zf\t,T,G,U)[s] = ^ Ps,^{t\s)Pu,G{s,X,T), (8) 

xec 
which depends on t via Ps,yL{t\s). 

As a consequence, the answer to a question like Preach{Zj^ ,t,T, G, U)[s] > p, p £ [0, 1], for Z^ 
depends on the initial time t: truth is time- dependent in Z\, . 



Example. We consider again the client-server example of Section 3.1 , and focus on two reachability 
probabilities for a single client: 

1. The probability of observing a time-out before being served for the first time within time T. 
This is a reachability problem with goal set G = {re} and unsafe set U = {rq,w}. 

2. The probability of observing a timeout within time T. This is a reachability problem with 
goal set G = {re} and unsafe set [/ = 0jj 



In Figures [5 (a) 5(b) , 6(a) and |6(b)] we can observe a comparison between the values computed for 



the limit ICTMC z and the exact ICTMC Z(^), for TV = 15 or iV = 150 (with a client-server ratio 
of 2:1), as a function of the time horizon T. As can be seen, the probability for z is in very good 
agreement with that of Z^^> (computed using a statistical approach, from a sample of 10000 traces) 
even for N relatively small. As far as running time is concerned, the fluid model checking is 100 
times faster for A'^ = 15, and 1000 times faster for N = 150, than the stochastic simulation. What 
is even more important is that the complexity of the fluid approach is independent of A^, hence 
its computational cost (on the order of 200 milliseconds for all cases considered here) can scale to 
much larger systems. Furthermore, another advantage of the fluid approach is that, by solving a 
set of differential equations, we are computing the reachability probability for each t E [0, T] (or 
better for any finite grid of points in [0, T]), while a method based on uniformisation (as in PRISM 
|16j ) has to deal with each time point separately. 

In Figures 5(c)[|5(d)||6(c) and 6(d), instead, we focus on the reachability probability for both 



problem 1 and 2 for T = 50 as a function of the initial time io £ [0,25]. The value for the fluid 
model is compared with the probability of Z^ obtained by simulating the full CTMC up to time 
to and then focussing attention on a specific client in state request and starting the computation 
of the reachability probability|j As we can see, the agreement is good also in this case. 

Finally, in Figures |5(e)[ |5(f)[ |6(e) and |6(f) , we compare the reachability probability for T 



100 (reachability problem 1) or T = 250 (reachability problem 2) of the ICTMC for different 
populations A^ and different proportions of clients (n) and servers {m), with the fluid limit. This 
data confirms that the agreement is good also for small populations for this model. 

""In fact, this is a first passage time problem. 

^This is done by using two indicator variables Xg and Xu that are set equal to one when a trajectory reaches a 
goal or an unsafe set, respectively. Then, we estimate the reachability probability by the sample mean of Xq at the 
desired time. 
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Figure 5: Client-Server model of Section [3. 1[ single client CTMC. First line: comparison of time- 
out before being served probability (point 1) for fluid and CTMC models as a function of time 
horizon T. Second line: comparison of time-out before being served probability (point 1) for 
fixed time horizon T = 50 and variable initial time Iq. Third line: time-out before being served 
probability (point 1) at time T = 250, and variable number of client and servers. 
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Figure 6: Client-Server model of Section 3.1, single client CTMC. First line: comparison of time- 



out probability (point 2) for fluid and CTMC models as a function of time horizon T. Second line: 
comparison of time-out probability (point 2) for fixed time horizon T = 50 and variable initial 
time to- Third line: comparison of time-out probability (point 2) at time T = 250, and variable 
number of client and servers. 
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5.2 Time- varying set reachability 

Now we turn our attention to the reachability problem for time- varying sets. First, we will focus 
on solving the problem for a generic ICTMC Z{t), considering then the limit behaviour of Z^ . 

In order to deal with the reachability problem for time varying sets, the main difficulty is that, 
at each time Tj in which the goal or the unsafe set changes, also the modified Markov chain that 
we need to consider to compute the reachability probability changes structure. This can have the 
effect of introducing a discontinuity in the probability matrix. 

In particular, if at time Tj a state s becomes a goal state, then the probability TTs-^^sit,Ti) 
suddenly needs to be added to the reachability probability from state si. Therefore, a change 
in the goal set at time Tj introduces a discontinuity in the reachability probability at time Tj. 
Similarly, if a state s was safe and then becomes unsafe, we have to discard the probability of 
trajectories that are in that state at time Tj, as those trajectories become suddenly unsafe. 

In the following, let G{t) and U{t) be the goal and unsafe sets, and assume that the set of time 
points in which G or U change value (at least in one state) is finite and equal to Ti <T2... < T^. 
This can be enforced by requiring that rate functions are piecewise analytic. Let Tq = t and 
Tk+i = t + T. 

In order to compute the reachability probability, we can exploit the semi-group property of 
the Markov process, stating that n(To,Tfc+i) = Y\,-^Qll(Ti,Ti+i). Then, we also need to deal 
appropriately with the discontinuity effects at each time Tj, mentioned above. We proceed in the 
following way: 

1. We double the state space, letting S = SUS, where a state s £ S represents state s when it 
is a goal state. Hence, in the probability matrix 11, '7rsi,s"2 is the probability of having reached 
S2 avoiding unsafe states, while S2 was a goal state. 

2. Consider a discontinuity time T and let ti G [Tj_i,Tj) and t2 G (Tj,Tj+i]. Define W{t) = 
S\{G{t)UU{t)). Then, for si G W{ti) and S2 G W{t2), the probability of being in S2 at time 
t2, given that we were in si at time ti and avoiding both unsafe and goal sets, can be written 
as 7fsi,52(ti,i2) = Y^seW{h)nW(t2)'^si,s{ti,Ti)7rs^s2iTi,t2). Hence, we have to appropriately 
restrict the summation set at time Tj, to account for changes in W. 

3. Consider again a discontinuity time Tj and let ti G [Tj_i,Tj) and ^2 G (Tj,Tj-|_i]. Suppose 
S2 G W{ti) and S2 G G{t2). Then, the probability of reaching the goal state S2 at time t2, 
given that at time ti we were in si, can be written as 

T^si,S2ih,Ti) + 2_^ Ti'si,s{tl,Ti)TTs,s2{Ti,t2)- 

sew(ti)nw{t2) 

The first term is needed because all safe trajectories that are in state S2 at time Tj sud- 
denly become trajectories satisfying the reachability problem, hence we have to add them to 
compute the reachability probability. 

All the previous remarks can be formally incorporated into the semi-group expansion of Il{t,t + 
T) by multiplying on the right each term n(Tj,Tj-|_i) by a suitable 0/1 matrix, depending only 
on the structural changes at time Tj+i. Let |5| = n and let Cvk(Tj) be the n x n matrix equal 
to 1 only on the diagonal elements corresponding to states Sj belonging to both W{T~) and 
W{Tj^) (i.e. states that are safe and not goals both before and after Ti), and equal to elsewhere. 
Furthermore, let Cg(^j) be the n x n matrix equal to 1 in the diagonal elements corresponding 
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to states Sj belonging to W{Tj^ ) n G(Tj^), and zero elsewhere. Finally, let C(^i) be the 2n x 2n 
matrix defined by: 



Consider now the following ICTMC Z on S, with rate matrix Q{t), where 

1. for si £ S and any S2 G S, gg^^sjlO = 0; 

2. for si W{t) and ah S2 G 5, g^si.saC*) = 

3. for si G W{t) and S2 G S* \ G(t), gsi.saC*) = Qsi,s2i't), while gsi,s2(*) = 0; 

4. for si G W{t) and S2 G G(t), qsi,s2it) = Qsusiit), while gsi,s2(i) = 0. 

In the previous chain, all unsafe and goal states are absorbing, while transitions leading from 
a safe state s to a goal state are readdressed to the copy s of s. States in S are absorbing, too. 

Now let n(ti, ^2) be the probability matrix associated with the ICTCM Q{t). Given the interval 
I = [t,t + T], we indicate with Ti, . . . , Tkj the ordered sequence of discontinuity points of goal and 
unsafe sets internal to /. Let 

T(t, t + T) = flit, Ti)c(ri)n(Ti, T2)C(r2) • • • c(Tfc,)n(T,,, t + t). (9) 

Then, we have that 

Ps{t) = Preach{Z,t,T,G,U)[s] = J^ T,,,-, (t, t + T) + l{s G G(t)}, (10) 

S1G5 

where the first term takes into account the probability of reaching a goal state starting from a 
non-goal state, while the second term is needed to properly account for states s G G{t), for which 
Ps{t) has to be equal to 1 (a formal proof can be given by induction on the number of discontinuity 
points). T(t, t + T) can be obtained by computing each n(Tj, Tj+i) solving the associated forward 
Kolmogorov equation and then multiplying those matrices and the appropriate ( ones, according 
to the definition of T. 

If we want to compute P{t) as a function of t, instead, we need a way to compute T(t, t + T) 
as a function of t. This can be done by observing that T depends on t only from the first and 
last factors in the multiplication. Defining r(Ti,Tfc) = (^(Ti)n(Ti, T2)C(T2) • • • n(T,fc_i,Tfc)C(Tfc), 
writing T{t,t + T) = Il{t,Ti)T{Ti,Tk)Il{Tk,t + T), differentiating with respect to t and applying 
the forward or backward equation for 11, we find the following differential equation for T: 

"^^^^dt^^^ " -Q(OT(t, t + T) + Tit, t + T)Q{t + T). (11) 

This equation holds until either t or t + T becomes equal to a discontinuity point. When this 
happens, the integration has to be stopped and restarted, recomputing T accordingly. 
Practically, to solve this problem we can proceed as follows: 

1. Given an interval [to) ^i] of interest for the initial time of the reachability, find all discontinuity 
points of the sets G and U contained in [to, ii + T], and let them be to = ^b < Ti < . . . < 
Tk < Tk+i = ti + T. Furthermore, let T/ = Tj + T for i = 0, . . . , /c, let pre{t) be the greatest 
Tj preceding t, and post{t) the smallest Tj following t. 
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2. Compute n(Tj,Tj-|_i) and Il{pre(T^),T[) for i < k, using the forward Kolmogorov equatfons 
Q Compute also each C{Ti)- 

3. Compute T(to,io + T) and integrate until time t = ioam{Ti,Tj^i — T}, where to + T G 

[Tj,Tj+i]. 

4. If t + T = Tj+i, multiply T on the right by ("(^j+i) and continue the integration. If t = Ti, 
then recompute T as n(ri,r2)r(r2, T,)li(rj,ri +r), where n(T,-,ri + r) = U{pre{T{),T[). 

5. Integrate piecewise using the previous rules until time ti. 

A more detailed algorithmic presentation of the procedure is given in Figure [7j Notice that, 
if the infinitesimal generator matrix Q{t) of Z is sufficiently well-behaved (for instance, Lipschitz 
continuous), then the function P{t) will be at least piecewise continuous, with a finite number of 
discontinuity points at instants Ti and T/. 

Remark 5.1. The precise behaviour of G and U functions at their discontinuity points (i.e. if 
they are left-continuous or right- continuous) is irrelevant for the computation of T: the set of 
trajectories of Z differing in those time points has probability 0. 

Remark 5.2. In the previous method, we need to integrate repeatedly a set An? differential equa- 
tions. However, most of these variables are redundant. In fact, we only need n^ variables for 
the probability transition matrix 11 on 5 and an additional n variables to store the reachability 
probability vector. The method presented above can be easily reconfigured to this restricted set 
of variables. 

Limit behaviour 

We consider now the limit behaviour of time- varying reachability probability for Z^ , proving 
that it converges (almost everywhere) to that of Zk- As in Section Hi we state this result in a more 
general form, assuming that also the goal and unsafe sets depend on A^, and converge robustly 
to some robust limit sets G and U. Furthermore, we require that G and U are compatible in 
the sense that they do not have a discontinuity at the same time for the same state s: Vs € 5, 
Disc{Gs) n Disc{Us) = 0. The following lemma, which is also the basic inductive step to prove 
convergence for CSL model checking formulae, relies on the functions involved being piecewise 
analytic. 

AN) 



Lemma 5.1. Let X^^' be a sequence of CTMC models, as defined in Section 



3.1 



and let Zj^ 



and Zk be defined from X^ ' as in Section 3.3, with piecewise analytic rates, in a compact interval 

[0,r'], for T' sufficiently large. 

Let G{t), U{t), t G [to, ti +T] be compatible and robust time-varying sets, and let G^^'{t), U^^'{t) 

be sequences of time-varying sets converging robustly to G and U , respectively. 

Furthermore, let P{t) = Preach{zk,t,T,G,U) and 

P(^)(t) = Preach{Zi''\t,T,G(^\U(^'>), t G [to,t,]. 

Finally, fix p £ [0, 1], MG {<, <, >, >}, and let Vpit) = I{P{t) ex p}, vi^\t) = I{P'^^\t) ex p}. 
Then 

1. For all but finitely many t G [to,ti]; -P^^'(t) — )• P{t), with uniform speed (i.e. independently 
oft). 



'^Notice, that, if Tj = pre{T[), then Yl{Tj,Tl) and n(rj,Tj+i) can be computed during the same numerical 
integration of the forward equation. 
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function reachability(Z, T, G, U, to, ti) 

Construct the CTMC on the modified state space S, according to the recipe in the text. 
Let to = To,Ti, . . . , Tfc, Tfc+i = ti+T be the time instants at which G or U has a discontinuity. 
for i = to /c do 

Compute n(Ti, Tj+i) and Il{pre{Ti + T), Ti +T) using the forward Kolmogorov equations 
end for 
Compute T{to,to + T) according to equation ^ and P{to) according to equation (10) 

t^to 

repeat 

Ta ^ pOSt{t) 

Tb ^ post{t + T) 
t^ min{ TajTfo -T] 



Compute T and P from t to t, according to ODE (11) and equation (10), with initial 
conditions T(t, t + T) (previously computed). 
if i + T = Tb then 

T{i,i+T) ^ T{i,i+T)C{Tb) 
else if t = Tq then 

T{i,i+ T) ^ U{Ta,post{Ta))r{post{Ta),pre{n))U{pre{n),n + T) 
end if 
t<^i 
until t > ti 
return [T(f , t + T),P{t)], t G [to, ii]. 
end function 

Figure 7: Algorithm for the computation of reachability probability P{t) for i G [to, ti] and time- 
varying goal and unsafe sets G{t) and U{t), with a finite number of discontinuities. Other input 
parameters are as in the text. 

2. For almost every p G [0, 1], Vp is robust and the sequence Vp converges robustly to Vp. 

Example. If we consider our running example, then it is easy to check that the rate functions 
defining the infinitesimal generator matrices of interest are piecewise analytic. In fact, even if the 
vector field of the fluid ODE is not analytic, due to the minimum function, the two functions 51 (x) 
and 92 (x) of which we take the minimum are analytic. Piecewise analyticity follows from the fact 
that the solutions of the associated ODE cross the surface fi'i(x) — 92 (x) = (where the minimum 
is not analytic) only a finite number of times. 

6 CSL Model Checking 

We turn now to consider the model checking of CSL formulae and the relationship between the 
truth of formulae for Zj, and Zk ■ 

Consider an until CSL formula (/? = Vtxip{ipi'U^^''^' (P2) , where (pi and ip2 are boolean combi- 
nations of atomic propositions. The major consequence of the time-inhomogeneity of z^ is that 
the truth value of (/? in a state s depends on the time t at which we evaluate such a formula. In 
particular, (p may be true in state s at time ti, but false at a different time ^2- Consequently, the 
set of states that satisfy a CSL formula (p can be time dependent, and this introduces an additional 
layer of complexity to the analysis of Zk- Indeed, this requires the computation of next-state and 
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reachability probabilities for time- varying sets. There is a similar problem with next formulae of 
the form ip = P^p(X' "'^'991), as the next-state probability also depends on the evaluation time t. 
Notice that we have the same issue about time-dependence also for model checking CSL formulae 
against Z^ (t). 

The method we put forward in the previous sections can cope with these issues, but in general 
may require a large computational effort (for until formulae, the solution of systems of ODE 
quadratic in the size of the state space of the ICTMC, and it depends on the number of discontinuity 
points of the sets U and G). However, in our setting we are interested in Zk, which is an abstract 
and approximate model of the behaviour of a single agent. Usually, a single agent has a very 
small state space, hence the given approaches to compute next-state probability and reachability 
of time- varying sets should be feasible in practice. 

An orthogonal issue is the asymptotic correctness of CSL model checking, when considering 
the sequence Zj^ and the limit Zk- As boolean operators pose no real problem, we only need to 
concentrate on next formulae ip = X' "' ''^992 and on until formulae ip = Vtxtpi^i'^^^ ' ^^2) , with 
time varying sets satisfying cpi and ip2- 

In particular, we can reduce this problem to the computation of the next-state probabilities 
P(^)(t) = Pnextizi^\t,Ta,Tb,G^^^) and Pit) = Pnextizk,t,Ta,Tb,G) (for next formulae) or to 

reachability probabilities pW(t) = Preach{zi'^\t,T,G^^\U^^^) and P{t) = Preach{zk,t,T,G,U) 

'or 
to 



(for until formulae), where G^ '{t) ([/' -'(t)) is the set of states satisfying (/?2 i~'^i ) at time t 



4.1 



and 



5.1 



Zj^ , while G and U are defined similarly for Zj. ^ Then, we may resort to Lemmas 
prove convergence of P^^>{t) to P{t) and of P^^{t) to P{t). 

However, in CSL model checking we are interested in truth values rather than in probabilities, 
and lifting the previous convergence to truth values is not so straightforward. Consider the path 
formula ipi\J'-'^''^'ip2. The problem is that we have to compute its probability P{t) (depending on 
the initial time t) for Zk and then solve the algebraic equation Ps{t) — p = for each state s, to 
identify for which time instants state s satisfies the formula. Now, the point is that, even in case 
p(^'(t) — )• P{t) uniformly, we are not guaranteed that P(^'(t) [xi p — )• P[t) ixi p. For instance, if 
P{t) = p, and IXI is <, then if P^ '{t) converges to P(t) from above, it never satisfies P^ '{t) ixi p 
for any A^, hence convergence of P^^\t) ixi p to P{t) ixi p does not hold. However, things can go 
wrong only when P{t) = p, and the main point of the convergence theorem is to prove that this 
happens sufficiently "rarely" not to impact on the computation of probabilities of a next or of an 
until formula in which tp is a sub-formula. 

In the following, we first outline an algorithm for CSL model checking of ICTMC, and then 
discuss convergence in more detail. Finally, at the end of the section, we will compare in more 
detail the CSL model checking problem for Z^ ' and {Zj^ \X(^'). 

6.1 Model Checking CSL for ICTMC 

The computation of next-state probabilities for time- varying target sets can be done by the method 
presented in Section |4j in particular the algorithm in Figure [3j 

The algorithm of Section |5.2| for computing reachability in the presence of piecewise constant 
goal and update sets, instead, is the core procedure to compute the probability of an until formula. 
In fact, consider the path formula ifiXJ^ "■'''' ip2. To compute its probability for initial time in 



*We can restrict our attention to until formulae with time between [0,r], as intervals [TajT;,] can be dealt with 
by essentially solving two reachability problems of this kind and combining their solution (or better, by computing 
two transient probabilities and then combining the two so obtained probabilities, see |13|). 
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[iO)*i]jjwe solve two reachability problems separately and then combine the results. 

The first reachability problem is for unsafe set Ui = [[-'</?il and empty goal set G{t+Ta) = 0. Let 
T^(t,t + Ta) be the probability matrix of this reachability problem. In order for the computation 
of the until probability to work, we must then discard the probability of being in an unsafe state, 



10 



essentially multiplying T^(t, t + Ta) by (^{t + Ta) on the right (see Section 5.2 ) 

The second reachability problem is for unsafe set U = [[-'V?i]] and goal set G = l(p2}, and is 
solved for initial time t S [to + Ta,ti + Ta], and time horizon T^ — Ta- Let T^(t + Ta,t + Tb) be 
the function computed by the algorithm in Section |5.2| for this second problem. Then, for each 
state s, safe at time t, we compute P{t) = T^{t, t + Ta)C^{t + Ta)T^(t + Ta, t + T}j)eg, where e^ 
is the vector equal to 1 for states s G 5 and zero elsewhere. Ps{t) contains the probability of the 
until formula in state s. Then, we can determine if state s at time t satisfies 'P^pl'/'iU^ "' *'c/92) by 
solving the inequality Ps{t) Mp. 

This provides an algorithm to approximately solve the CSL model checking for ICTMC recur- 
sively on the structure of the formula, provided that the number of discontinuities of sets satisfying 
a formula is finite and that we are able to find all the zeros of the computed probability functions, 
to construct the proper time-dependent satisfiability sets (or approximations thereof). The full 
procedure is sketched in Figure |8J 

Below we will consider this algorithm in more detail, focussing particularly on correctness and 
termination. In this consideration we will make the following assumption about the numerical 
algorithms that it uses. 

Assumption 1. There are interval arithmetic routines that can compute bounding sets for the rate 
functions of z^ and Zj^ , in such a way that the approximation error can be made arbitrary small. 
We call such functions interval computable. 

Notice that this assumption is not very restrictive. It applies to all the standard functions, 
and also to solutions of ODEs of functions which satisfy it, to derivatives of these functions and 
to their integrals |45t I46j . In particular, if the rate functions are interval computable, then so will 
be all the probabilities computed by solving reachability problems. 

The approach presented above relies, in addition on the solution of ODEs, also on two other 
key numerical operations: given a computable real number p, determine if p is zero and and given 
an analytic function /, find all the zeros of such a function (or better an interval approximation 
of these zeros of arbitrary accuracy). However, it is not clear if these two operations can be 



carried out effectively for any input that we can generate, see Remark 6.2 for further comments. 
Therefore, we need some further assumptions. Instead of restricting the class of functions (which 
seems a difficult problem since we have to consider the solution of differential equations), we will 
follow the approach of [17], introducing a notion of robust CSL formula and proving decidability 
for this subset of formulae. This will not solve the decidability problem in theory, but makes it 
quasi-decidable |37], which may be enough in practice. As we will see, the set of CSL formulae 



which is not robust has measure zero (see Theorem 6.2). 

In order to introduce the concept of robust CSL formula, consider a CSL formula (p and let 
pi, . . . ,Pk be the constants appearing in the Vt^jp operators of next and until sub- formulae of c/3. 
We will treat f = (p{pi, ■ ■ ■ ,Pk) as a function of those pi, . . . ,pf^. Furthermore, we will call the 
next or until sub-formulae of ip top next sub-formulae or top until sub-formulae if they are not sub- 



^The appropriate value of to and ii are to be deduced from (pi, if>2 and the superformula of the until, in a standard 
way [44] 

^°In fact, this reachability problem can be solved in a simpler way: it just requires trajectories not to enter an 
unsafe state, and then collects the probability to be in a safe state at the time t + T. In particular, we can get rid 
of the copy 5 of the state space, and define a simplified T function using (^w matrices instead of ^ ones. 
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function CSL_MC(Z, if, to, ti) > Computes Vs{t) = I{s,t \= (p} for s £ S and t £ [to,ti]. 

if (p = p then 

Vs{t) ^ l{p G L{s)], s G S 
else if (/3 = ^Lfi then 

Vi ^CSLJMC{Z,(pi,to,ti) 

V{t)^l-Vi{t) 
else if ip = (fi f\ (p2 then 

Vi ^CSL_MC(Z,(/7i,to,ti) 

V2 ^ CSh_MC{Z,^2MM) 

y(t)^mm{Vi(t),y2(t)} 
else if v? = Pc^plXl^^'^^Vi) then 

Vi ^CSL_MC(Z,(/7i,to,ti) 

P <r- next-state-probability(Z, Vi, Tq, T5, io, ii) 

y(t)^i{p(t)Mp} 

else if (/9 = Ptxip(v'iU["^"'"^''V2) then 
Vi ^ CSL_MC(Z, -v^i, to, ti + Tb) 
V2 ^ CSL_MC(Z, (^2, to, h + Tb) 
T^ ^ reachability(Z, Tfl, 0, ¥l, to, ^i + Ta)[l] > Returns T component of 

REACHABILITY. 

T2 ^ reachability(Z, n - Ta, V2, Vi,to + Ta, h + Tb)[l] 
P{t) = T\t, t + Ta)C\t + r,)T2(t + T„, t + rfe)e5- 
V{t)^l{P{t)ixp} 
end if 
return F 
end function 

Figure 8: Core algorithm for solving the CSL model checking problem, by computing the satisfia- 
bility of a CSL-formula 99 as a function of the time t G [to, ti] at which it is evaluated. The truth 
value of if is then the value it has in to, which is usually 0. 

formulae of other next or until formulae. The other next or until formulae will be called dependent. 
Finally, given two robust time- varying sets Vi and V2, we recall that Vi and V2 are compatible if 
they do not have discontinuities for the same state s happening at the same time instant t. 

Definition 6.1. A CSL formula 99 = y^{p), p G [0, 1] is robust if and only if 

1. there is an open neighbourhood H^ of p in [0, 1]^ such that for each pi G W, 

■5,0 1= (p{p) <^ S,0 1= (/9(pi). 

2. The time-varying sets of any dependent next or until sub-formula of ip holds are robust. 

3. The time-varying sets of sub-formulae of if that are part of an until formula or of a conjunc- 
tion/disjunction are compatible among them. 

We now prove the following theorem, which states that the CSL model checking algorithm we 
put forward works at least for robust formulae: 

Theorem 6.1. The CSL model checking for ICTMC, for piecewise analytic interval computable 
rate functions, is decidable for a robust CSL formula ^(j>i, . . . ,Pk)- 
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The following corollary is a straightforward consequences of the proof of the previous theorem: 

Corollary 6.1. The algorithm for CSL model checking presented in this section is correct for 
robust CSL formulae. I 

We turn now to characterise the set of robust formulae from a topological and measure-theoretic 
point on view. We have the following 

Theorem 6.2. Given a CSL formula ^{p), with p € [0, 1] , then the set {p | ^{p) is robust} is 



relatively oper in [0, 1] and has Lebesgue measure 1. 



The openness of the set of robust thresholds for a formula allows us to prove the following 
corollary about quasi- decidability. In this paper, we consider a notion of quasi-decidability, which 
is slightly different than the one defined in [33 • In fact, we take advantage of the fact that out 
input values belong to a compact subset K C M", for which a standard notion of measure exists. 

Definition 6.2. A problem with inputs in a compact subset K C M" of Lebesgue measure m{K) > 
0, is quasi- decidable if there is an algorithm that solves it correctly for an open subset U G K, 
with fie{U)/fie{K) = 1. 



Combining Theorems 6.1 and|6.2[ we obtain the following: 



Corollary 6.2. The CSL model checking for ICTMC, for piecewise analytic interval computable 
rate functions, is quasi- decidable for any formula ip{p). ■ 

Remark 6.1. The notions of robustness and quasi-decidability have a practical side. First, the 
openness property of the set of robust thresholds for a formula ip{p) guarantees that if we "perturb" 
a formula (by varying the set p of threshold constants of the path probability operators), then 
the formula remains robust. Furthermore, by the definition of robustness, also its truth value 
remains the same (as the notion of quasi-decidability of [42] requires). This explains the use of the 
terminology "robust" . 

Secondly, the characterisation of the set R of robust thresholds for a formula ip provided in 



Theorem 6.2, implies that if we choose thresholds at "random", we are likely to select a robust 
formula. In fact, consider the grid of rational numbers with - in [0,1], i.e. GRn = {jj; \ it^ < 
n, rra, n G N}, and take the Cartesian product Gi?^ C [0, 1]*^. Let /U„ be the uniform distribution in 
Gi?^, then ij,n — )■ M) the uniform distribution on [0, l]'^ (which coincides with the Lebesgue measure 
on Borel sets). Now, as R is open and has Lebesgue measure 1, then ^{R) = 1 and fi{dR) = 0, 
hence Ris a continuity set for fi. Therefore, Hn{R) — ^ A'(-R) = 1 by the Portmanteau theorem 



This means that, fixing e > 0, if we choose the thresholds of the until sub-formulas from the set 



GR^, for n large enough, the probability of choosing a bad set of thresholds, for which the formula 
is not robust and the CSL model checking algorithm may not terminate, will be less than e. 

Remark 6.2. The semi-decidability result presented here is in contrast with the decidability result 
of model checking for time- homogeneous CTMC. However, in that case the result follows because 
^,(0) has a special form allowing the application of Lindeman-Weierstass theorem for transcen- 
dental numbers (together with zero testing procedures for algebraic numbers) |49j . This, in turn, 
is a consequence of having constant (rational) rates. In our case, instead, rates are piecewise 

^^ A set U G V is relatively open in 1/ C W, where W is a. topological space, if it is open in the subspace topology, 
i.e. if there exists an open subset (7i C W such that U = V nUi. 
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analytic functions, and we cannot rely on the method of [|49] anymore. In fact, in the algorithm 
for computing the probability, there are two numerical operations that are potential sources of 
undecidability: 

1. Given a number p, which is the analytic image of a rational, decide if it is zero. This is a 
classical problem whose decidability is not known, even restricting to expressions made up 
by polynomials and exponentials only |50tl51j. Indeed, its decidability is connected with the 
truth of the Schanuel conjecture |50ll51j . which is in turn connected with decidability of the 
theory of reals extended by the exponential. However, even in case the Schanuel conjecture 
holds, it is not clear if the zero problem will be decidable for any analytic function. 

2. Detecting the zeros of an analytic function with arbitrary precision. In this case the problem 
is caused by non-simple zeros, i.e. points in which the function and some of its derivatives 



are zero. The method sketched in a footnote of the proof of theorem 6.1 does not work, as 
it relies on the fact that we can bound the derivative away from zero on null points of the 
function. Furthermore, in the presence of non-simple zeros, detecting if a compact interval 
is bounded away from zero is semi-decidable (the decision procedure fails if the interval 
contains a non-simple zero). Whether there is a decidable algorithm for this problem is not 
known to the authors (even assuming the Schanuel conjecture is true). It may be possible, 
however, to find algorithms for some subclass of analytic functions large enough for practical 
purposes. For instance, if we know a lower bound on the radius of convergence of power 
series in each analytic point, we can effectively extend the real analytic function to an open 
ball in the complex plane, and then use methods developed for complex analytic functions 
|52j which can effectively compute the number of zeros in any sufficiently simple open set, 
by integrating a function on its boundary with interval arithmetic routines \53\ [52] . 

Our conjecture is that the model checking problem for time-inhomogeneous CTMC is not 
decidable in general, although it may be decidable for some restricted subclass of rate functions if 
the Schaunel conjecture is true. Further investigations on this issue are required. 

Finding an upper bound on the complexity of the approximation algorithm, when it converges, 
requires us to find an upper bound on the number of zeros of the solution of a linear differential 
equation with piecewise analytic rates. This is a non trivial problem. However, we can rely on a 
result for linear systems with bounded analytic rate functions [54], which gives an upper bound ^ 
on the number of zeros, expressible as an elementary function of the upper bound on coefficients 
of the ODE. For piecewise analytic rate functions, simply multiply this bound for the number of 
analytic pieces. The number of analytic pieces is Kq + K, where Kg comes from the piecewise 
analytic nature of rate functions, and K from the number of structural changes of [[~"/^iI| and [[c/52]] 
sets. Hence, the number of zeros of -P(t) — p can be bounded by {K + Kq)'^. By induction, if h 



is the degree of a formula (p and hu is the number of nested next or until subformulae, ^^ then the 
total complexity is bounded above by C{n,Tmax:£){'^^'^^'^KQ), where the constant C{n,Tmax,£) 
hides the cost of integrating ODEs and finding roots in each analytic piece. It is proportional to 
n^ (matrix multiplication), time T^ax for which the ODEs have to be solved, and the precision e 



of root finding and numerical integratior ^^ 



However, this is a theoretical upper bound, and we will not expect such a complexity in practice. 



^■^Notice that, for next or until formulae not containing any other next or until subformula, K = 
^^The precision e depends on the specific analytic functions considered. However, we can imagir 
which takes an e as input and does not provide an answer if the precision is not small enough. 
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6.2 Convergence for CSL formulae 

We are now ready to state a convergence result for CSL model checking. Also in this case, we will 



restrict our attention to robust CSL formulae. This is reasonable, as we want to use Lemmas 4.1 



and 5.1, which require robustness of time- varying sets. 



Theorem 6.3. Let X^^' be a sequence of CT MC models, as defined in Section 
and Zfc be defined from X^'^' as in Section 



3.3 



3.1 



and let Zj^ 



Assume that Zj^ , Zj. have piecewise analytic infinitesimal generator matrices. 

Let (p{pi, . . . ,pk) be a robust CSL formula. Then, there exists an Nq such that, for N > Nq and 

each s £ S 

s,Ot= (iv) if 4^ s,0^zk f- 

Corollary 6.3. Given a CSL formula ^{p), with p € [0,1]^^', then the subset of [0,1]^^ in which 
convergence holds has Lebesgue measure 1 and is open in [0, 1] . ■ 

The previous theorem shows that the results that we obtain abstracting a single agent in a 
population of size N with the fluid approximation is consistent. However, the theorem excludes the 
sets of constants p for which the formula is not robust. Interestingly, this is the same condition re- 
quired for decidability of the model checking problem for ICTMC, a fact that shows how these two 
aspects are intimately connected. Notice that, contrary to decidability, this limitation is unavoid- 
able and is present also in the case of sequences of processes converging to a time-homogeneous 
CTMC. In this case, in fact, the next-state and reachability probabilities are constant with respect 
to the initial time, and their value p (in the limit model) can cause convergence of truth values to 
fail. 

However, notice that the constants p appearing in a formula that can make convergence fail 
depend only on the limit CTMC z^, hence we can detect potentially dangerous situations while 
solving the CSL model checking for the limit process (in these cases the model checking algorithm 
may fail to provide an answer). 

Remark 6.3. In this paper, we are considering only time bounded operators. This limitation is 



a consequence of the very nature of the approximation theorem 3.2 which holds only for a finite 
time horizon. However, there are situations in which we can extend the validity of the theorem to 
the whole time domain, but this extension depends on properties of the phase space of the fluid 
ODE [SSI IMl Ell- 
in those cases, we can prove convergence of the steady state behaviour of Zf^ to that of Zk, 
hence we can incorporate also operators dealing with steady state properties. 

In order to deal with time unbounded operators, instead, convergence to steady state is not 
enough. We also need to ensure that the equation P{t) — p has a finite number of zeros on the 
whole positive time axis. Piecewise analyticity is not sufficient in this case (think about sine and 
cosine), and stronger conditions have to be required. However, for periodic functions, we may 
reason similarly to [26J, if we can prove that periodicity of rate functions implies periodicity in the 
reachability probabilities as a function of initial time. 

Example. Going back to the running example, consider the until path formula trueU^^'^^' timeout, 
where timeout is true only in state re. Its probability, as a function of the initial time, is shown 
in Figures 9(a), |9(b)[ and |9(c) for the states rq, w, and t, respectively. In the same figures. 



we also show the time- dependent truth of the CSL formula V<^Q,iQT{trueU^^'^^^ timeout), which is 
obtained by solving the inequality Ps{t) < 0.167, where Ps{t) is one of the previous time-dependent 
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probability functions. In this case, we can observe that for time to ^ [0,100], there is only one 
solution, as the probability is monotone. This depends on the solution of the fluid equations. In 
this case, in fact, they converge to a steady state, hence we do expect that also the time dependent 
truth value of CSL until formulae stabilises (when the fluid ODE are close to steady state, the 
rates of the ICTMC are practically constant). This suggests that in many practical cases, the 
number of changes of truth value of until formulae will be very small, as in the running example. 
Notice that in the case of the running example, if we had chosen a threshold bigger, say, than 0.25, 
then the time-dependent truth formulae would have been a constant function. 
In Figure [9j instead, we show the probability of the path formula 



tTOeUl°'^l(P<o.i67(iraeC/[°'^°UimeoMt)), 



as a function of the time horizon T. In the plot, it is evident how this probability has discontinuities 
in those time instants in which the truth values function of its until sub-formula change. These 
discontinuities differentiate the model checking of ICTMC from that for time-homogeneous CTMC. 



Pn^[F<=iO timeout] ^tO varying 



Pn=?[F<=*0 timeout] —tO varying 





(a) iraeU[°''^°l timeout- 



rq 



(b) true'U^°'^°Uimeout - w 



Pn^[F<=40 timeout] —tD varying 



(c) irueU'"'^") timeout - t 



(d) iraeU[°'^l(-p<o.i67(i™eU[°'^°l timeout)) 



Figure 9: Figures 9(a) , 9(b) and 9(c) Probability of the formula true Ul'^'^'^l timeout, for varying 
initial time, and different initial states {rq, w, and t respectively). The dotted line shows the 
time varying truth function for the CSL formula V<^o,ie7{true\J'^^'^^' timeout)), which is obtained 
by finding the zeros of the initial-time dependent probability. Figure |9(d)[ Probability of the until 
path formula trueU['^''^](P<o.i67(^'"weU['''^'^] timeout)), as a function of time bound T. Vertical 
dotted lines show the discontinuity points time-dependent truth of the until sub-formula. 
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6.3 Comparison of CSL model checking for Z^. and {Z^ ,X*^^^) 

In this paper we have considered two possible descriptions of a single agent at a fix ed popu 
level N, i.e. Z^. (t) and {Z^ {t),'K^'^'{t)). From the discussion in Sections 3.3 4, and Sl 



ation 
we 



already know that, while {Zj^ {t),'K^^'{t)) is a CTMC with finite (but extremely large) state 
space, Zj^ (t) has a much smaller state space but it is not a Markov process. Furthermore, 



W/ 



its behaviour is time dependent. The non-Markovian nature of Zj;. (t) has consequences for its 



reachability probability (see Section 5.2), meaning that its value is dependent on the initial time at 



which we compute it. This implies that the satisfiability of a CSL formula (with the truth value of 
atomic propositions depending only on S) for Zj^ (t) can depend on the time at which we evaluate 
it. Hence we need to consider time-dependent sets to compute the probabilities of next or until path 
formulae. But time-dependent sets can introduce discontinuities in such probabilities, as discussed 
On the other hand, {Zj^ (t),X( ■'(t)) is a time-homogeneous CTMC, hence its 



5.2 



in Section 

next-state and reachability probabilities do not depend on time and no time-dependent notion of 
satisfaction has to be considered in this case. In particular, when considering (Z^ (t), X(^-'(t)), 
its reachability probability is always a continuous function. This implies that the truth value of 
a CSL formula containing nested next or until sub-formulae, can be different if we consider its 
satisfiability with respect to zl^\t) or (Z^^^(t), X(^)(t)). 

However, despite this discrepancy for finite A^, we will prove that the satisfiability for Zj^ (t) 
and {Zj^ (t), X''^'(t)) is asymptotically the same, at least if we restrict to robust CSL formulae. In 
order to show this, we will combine the convergence results of the previous sections with additional 
results relative to {zf\t),±^^\t)) and (z(t),x(t)). 
Example. If we observe Figures 6(c)| and |6(d) we can easily convince ourselves that the reachability 



probability for Zj^ in the running example for the formula ipi = trueXJ'^'^^' timeout depends on 
the initial time. Hence it gives rise to a time-dependent set for the satisfiability of the formula 
ip2 = 'P<o. 167(951)- This implies that for Zj^ , the probability of the formula if = trueXJ'^' ' (/?2 will 
have discontinuities as a function of T, similarly to the case for z^- However, if we compute the 
reachability probability for ip in (Zj^ ,X'^''), in a state s,xo, this will be a continuous function 
of T, hence the two probabilities are different. 

We will now prove the convergence of the standard CSL model checking for Y^ •'(t) = (Z^ (t), 
-^ (^)) ™ state s,xo, to the equivalent CSL model checking procedure for y(t) = (z(t),x(i)). 
This procedure requires us to compute, given a next formula ip = Xl^'"^'']^?! or an until formula 
If = ipiXJ'- "^^ '''ip2, its probability P(s,x) starting from time 0, in each point (s,x) of the state 
space S X E oi y(t), and then solve the inequality P{s, x) ixi p, to determine the truth of V^p{ip) 
in (s,x). This defines a subset of 5 x £^ where 'P^p{ip) is true. 

The intuition behind the proof is that the truth value of an until formula in a state (s,xo) for 
y(t) does not depend on the whole state space S x E, but only on the points of E intersected by 
the solution of the fluid ODE starting in xq, i.e. on 5 x $([0,T],xo), where $(t,xo) is the flow of 



the differential equation ^^ Furthermore, the convergence of X^ ■* (t) to x(t) allows us to restrict 
the attention to an arbitrary small neighbourhood of <l>([0,T],xo), in order to solve the model 
checking problem for Y(^-*(t), for N large enough. 

In the following, we need some additional concepts and definitions. 

Consider the domain P'^-* C E of X*-^-*. With each point x G i?, we associate a point 
z/*-^)(x) G V^^\ such that ||x — z^'^'*(x)|| < jj. The existence of such a point is guaranteed by 



*The solution of the fluid ODE at time t starting in xq at time 
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the definition of E. Now, we further assume that, given a point (s,x) G E, the initial state 
Y'^(O) is {sjU^ '{'k)), so that Y'^(O) converges to (s,x) uniformly in space. This choice of 
Y(^)(0) guarantees uniform bounds in space for Kurtz theorem and the fast simulation theorem, 
for convergence in probability}^ 

Now, consider the fluid limit differential equation, and let <I>(t, xq) be its flow. We assume that 
$(t,xo) is a piecewise analytic function with respect to t and x. The T,£-flow tube for xq is the 
set £^0 C E, defined by £^0 = *^([0, T], Bs{xo)), i.e. the set of all trajectories up to time T starting 
in a ball of radius e centred in xq. Now, consider a T, e-flow tube Eq for xq. For any x £ Eq, let 
T+ = T+(£'o) = sup{t I ^{[0, t],x) G Eq} be the time at which the trajectory starting in x leaves 
Eq. Furthermore, let T^ = T^[Eq) = inf{i | <I>([t, 0],x) G £"0} be the time at which the trajectory 
starting in x enters Eq. 

A subset L> C 5 X £0 is a d-set for Eq if and only if, (i) D is closed (in S x £9), (ii) D is 



the union of a finite number of smooth manifolde of dimension n — 1 or less, and (iii) for each 



X G £0, it holds that {s} x <l>([T~(£o),T+(£o)],x) n D contains at most k points in each state s. 
In other words, a d-set is a union of piecewise analytic manifolds that intersects each trajectory in 
at most k points. It can be easily checked that each d-set has (Lebesgue) measure zeroFj 

We also introduce a notion of robust subset of 5 x Eq, for a T, e-flow tube £0 in xq. Consider a 
subset V C S X Eq. We say that V is robust in S x Eq if and only if, (i) its boundary dV is a d-set 
in 5 X £0, and (ii) for each (s,x) G 5 x Eq, the time-varying set l^[s](t) = l{(s, <I>(t,x)) G V}, 



T^ < t < T^, is robust in the sense of Definition 4.2 (notice that it contains at most k < 00 
discontinuity points, where k does not depend on x, as dV is a d-set). We sometimes denote 
dV by Disc{V). We also say that two robust subsets Vi and V2 of 5 x Eq are compatible if 
dVi n dV2 = 0. 



Similarly to Section 5^, we say that a sequence of sets V^'^' C S x Eq converges robustly to a 
robust set V '^ S x Eq, with Eq a T, e-flow tube in xq, if and only if, for each open neighbourhood 
U of Disc{V), there is iVo > such that, ViV > iVo and ah (s,x) G (5 x £0) \ U, (s,x) G F^^^ if 
and only if (s, x) G V. 



We are now ready to state the following lemmas, which are space- versions of Lemmas 4.1 and 



5.1 on time-varying sets, and are the key to the induction step of Lemma 6.4 



Lemma 6.1. Let Eq C E be a T,eQ-flow tube for xq. Let G be a robust subset of S x Eq, and 
G^^' be a sequence of subsets of S x Eq that converge robustly to G. 

Let P{s,x) = Pnext{'y,s,x,Ta,Tf,,G) be the probability that the first jump ofy{t) is into a state 
in G and happens at a time t G [T^, Tf,], given that y started at time t = in state (s, x) £ S x Eq, 
and let P^^'{s,x) = P^^J^{Y^^',s,u^^'{x),Ta,Th,G^^') be defined similarly, withG andx replaced 
by G^^' and z^(^'(x), respectively. Furthermore, define V = {(s,x) | P(s,x) cxi p} and 

y(^) = {(s,x) I P^^'{s,x) txip}. Then there exists ei > such that, in £1, the (T — Tf,),£i-fiow 
tube for xq: 

1. P^ '{s,x) —7- P{s,x) for all x G £1, uniformly in (s,x). 

2. IfVx:Q{t), t G [^x^, (£i),7'x'o(-^i)]? ^-^ '^ robust time-varying set, then V is robust in £1 and 
V^ ' converges robustly to V. 



^^The speed of convergence to the fluid limit depends on the initial conditions only through ||X'-^'(0) — x(0)||; the 
choice of u^^' (x) guarantees the uniform convergence of this quantity with respect to x. 

^''An smooth manifold is the zero set of a sufficiently smooth function, in this paper at least having continuous 
first-order derivatives. 

^'^ Any set of topological dimension n ~ 1 or less has Lebesgue measure zero in R". 
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Lemma 6.2. Let Eq C E be a T,eQ-flow tube for xq. Let U and G two robust and compatible 
subsets of S X Eq, and U^ ' , G^ ' be sequences of subsets of S x Eq that converge robustly to U 
and G, respectively. 

Let P{s, x) = Preachiy, s, X, Ti, T2, U, G) be the probability that y(t) reaches a state in G within 
time [Tq, Tf,], avoiding any unsafe state in U, given that y started at time t = in state (s, x) G 5 x 
Eo, andletP(^\s,x) = PJ:^J^^{Y'^^\s,u'^^'){x),Ta,n,U'^^\G^^^) be defined similarly, withG, U, 
X replaced by G^ \ U^ ' , andu^ '{yi), respectively. Furthermore, define V = {(s,x) | P(s,x) ixi p} 
and y(^) = {(s,x) | P^^'{s,:x.) txip}. Then there exists ei > such that, in Ei, the (T — Tf,),ei- 
flow tube for xq ." 

1. P^ '{s,:x.) —7- P{s,'x.) for all Ji. £ Ei\ D, where D is a d-set, uniformly in (s,x). 

-2- // V"xo(i), t G [r~^(£'i),T+ (£'i)], is a robust time-varying set, then V is robust in Ei and 
V^ ' converges robustly to V. 

The previous lemmas are the key arguments used in the structural induction to prove the 
following result. 



Lemma 6.3. Let X^^' be a sequence of CTM C models, as defined in Section 
and Zk be defined from X^ ' as in Section 



3.3 



3.1 



and let Z 



(N) 



Assume that there is a flow tube Eq of 'Kq such that all trajectories in Eq are piecewise analytic. 
Let if = 9j(p) be a robust CSL formula for the trajectory <^(t,xo). Then, there is an Nq such that, 
for all N >Nq, 

S,X0 \=y ip <^ S,V^^\yiQ) \=y;{N) ^■ 

We now turn to consider the relationship between the model checking problem of a CSL formula 
if for Zk{t) and the model checking problem for the same formula with respect to y{t). In this 
case, it is easy to see that a formula is true for Zk{t) if and only if it is true for y(t). In fact, in this 
process the truth value of a formula in state (s, xq) depends only on the trajectory <&(t, xq) starting 
in Xq. Furthermore, if we fix a time i and consider the point Xf = <I>(t, xq), then the process Zk{t), 
defined with respect to the trajectory $(t,Xf) starting in point Xf at time zero, equals the process 
Zk{t + t), starting in xq at time zero, due to the semi-group property of the flow $(•, •). Hence, any 
reachability probability for z^ with respect to the initial time t equals the reachability probability 
for Zk at time 0: We can always turn a time-dependent reachability problem into a more classical 
space-dependent one. From the previous discussion, the following lemma follows: 



3.1 



and let Z 



(N) 



Lemma 6.4. Let X^ ' be a sequence of CTM C models, as defined in Section 

and Zk be defined from X^'^' as in Section 3.3. 

Let if = 93(p) be a robust CSL formula for the piecewise analytic trajectory $(t, xq), and let Zk be 

the ICTMC defined on S with respect to trajectory $(i,xo). Then, 



Using the previous lemmas, we can easily show the following theorem. 



Theorem 6.4. Let X^^' be a sequence of CTMC models, as defined in Section 
there is a flow tube Eq of xq such that all trajectories in Eq are piecewise analytic. 



3.1. Assume that 
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Let Lp = v^Ip) be a robust CSL formula for t he tr ajectory <I>(t,xo), let Zj^ (t) and Zk{t) be the 



stochastic processes on S defined as in Section 3.3. and let y{t) and Y^ '{t) be defined as in this 



section. Then, there is an Nq such that, for all N > Nq, 

S \= (N) ^ <^ S,Z^(^)(xo) \=Y(N) (p. 



k 



Proof. There exists an Nq, such that, for all N > Nq, 

^(N) if <^ S\=zk -^ S,Xo|=yC^ <^ S,iy'-^'>{:S.Q) \=Y(N) if, 



^k 



where the first equivalence follows from Theorem 6.3, the second equivalence from Lemma 6.4, and 



the third equivalence from Lemma 6.3, while Nq can be chosen as the largest one between that of 
Theorem 16.31 and that of Lemma 16.41 ■ 



Inspecting the proof of the previous theorem, the following corollary is straightforward. 

Corollary 6.4. Let ip = (p{p) be a robust CSL formula for the trajectory <l'(t,xo). Then, there is 
an Nq such that, for all N > Nq, 

S,U^^\-Kq) |=Y(iV) if <^ S \=z^ if. 



7 Conclusions 

In this paper we exploited a corollary of fluid limit theorems to approximate properties of the 
behaviour of single agents in 1 arge population models. In particular, we focussed on reachability 
and stochastic model checking of CSL formulae. The method proposed requires us to model check 
a time-inhomogeneous CTMC of size equal to the number of internal states of the agent (which is 
usually rather small). Hence, it gives a large improvement in terms of computational efficiency. 

We then focussed on the reachability problem for ICTMC, both in the case of time-constant and 
time-varying sets. We provided algorithms to tackle both cases, and we also proved convergence 
of the reachability probabilities computed for the single agent in a finite population of size A'^ to 
those of the limit fiuid CTMC. Finally, we focussed on model checking CSL formulae for ICTMC 
proposing an algorithm that works for a subset of CSL including the time bounded next operator 
and the time bounded until operator. We also showed decidability and convergence results for 
robust formulae, proving that the set of non-robust formulae has measure zero. 

There are many issues that we wish to tackle in the future. First, we would like to better 
understand the quality of convergence. This can be accomplished by trying to derive theoretical 
error bounds (which may be too loose to be of practical interest) and by running many experiments 
to identify situations in which the approximation performs well (in terms of both classes of formulae 
and model structure). In addition, we would like to provide a working implementation of the model 
checking algorithm for ICTMC, studying its computational cost empirically (and exploring how 
easy it is in practice to find a non computable instance). Furthermore, we want to investigate the 
connections between single agent properties and system level properties. We believe this approach 
can become a powerful tool to investigate the relationship between microscopic and macroscopic 
characterisations of systems, and to understand their emergent behaviour. 
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As far as CSL model checking for ICTMC is concerned, we aim to extend it to include time 
unbounded and steady state operators, at least for those subsets of rate functions in which the 
algorithm can be shown to be decidable. We also need to consider rewards, at least for a finite 
time horizon (here we expect their inclusion to be relatively straightforward). Then, we would like 
to show convergence results also for this larger subset of CSL, under the hypothesis required for 
steady state convergence of the fluid approximation. 

Another line of investigation would be to consider different temporal logics, such as MTL. 
For this logic, asymptotic correctness is relatively easy to prove, along the lines of Proposition 



5.1 What is more difficult is to find an effective algorithm to model check MTL properties for 
ICTMC. One possibility may be to combine the approaches of |3 H I26 | [27]. and exploit algorithms 
and techniques to compute reachability of PDMP 
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A Proofs 

In this appendix, we present the proofs of propositions, lemmas, and theorems of the paper. 
A.l Next-State Probability 



Proposition (4.1). Let /:/—)• R 6e a piecewise analytic function, with / C R a compact interval. 
Let Ef = {x ^M.\ /i£(/~^({x}) = 0} be the set of all values x such that f is not locally constantly 
equal to x, where fii is the Lebesgue measure. Furthermore, let Z^ = f~^i{x}) be the set of 
solutions of f{t) = X and let DZf = {x G R | Vt G Zx, fit) ^ 0}. Then 

1. \/x G Ef, Zx is finite. 

2. Hi{Ef n DZf) = 1 

Proof. Point 1 follows from basic properties of the piecewise analytic function (/ — x): in any 
analytic piece, either the function is constantly equal to zero, or it has only a finite number of 
zeros. Point 2, instead, follows from the fact that the derivative f'{t) of t is piecewise analytic, 
hence has only a finite number of zeros (in the analytic pieces in which / is not constant). ■ 



Proposition (4.2). Let V^'^'lt) be a sequence of time varying sets converging robustly to a robust 



set V{t), t G /. Let L'Jf ^ = {t \ V^^\t) / V{t)}. Then Hi{d\^^) -^ 0, where m is the Lebesgue 



measure on 



Proof. A straightforward consequence of the definition of robust convergence is that, for each open 
neighbourhood U of DisciV), there exists an A'^o such that, for all N > Nq, V^^\t) = V{t) 
for t G I \ U. Now, as V is robust, then \Disc{V)\ = m < cxd. Fix e > and define 
Us = UteDjsc(v) -^(^' ^)' where B{t, e) is the open ball centred in i of radius e. Then ^i{Ue) < 2me. 
Now, fix Ek -^ 0. For each /c, there is an Nk such that, for all N > Nk, V^'^\t) = V{t) for t G I\Ue^, 



and therefore Dy C U^^. ■ 



Lemma (4.1). Let X^^' be a sequence of CTMC models, as defined in Section 



3.1 



and let Z^ 



and Zk be defined from X^ > as in Section 3.3, with piecewise real analytic rates, in a compact 



interval [0, T'], for T' >ti + Tb. 

Let G{t), t G [to,ti+Th] be a robust time-varying set, and let G^'^'it) be a sequence of time-varying 

sets converging robustly to G. 

Furthermore, let P{t) = Pnext{zk,t,Ta,Th,G) and P^^>{t) = 

Pnext{Zf. ,t,Ta,Tb,G), t G[tQ,ti]. 

Finally, fix p G [0, 1], MG {<, <, >, >}, and let Vp{t) = /{P(t) ex p}, V^^\t) = /{pW(t) cxi p}. 
Then 

1. P'^^^t) -^ P{t), uniformly in t G [to,ti]. 

2. For almost every p £ [0, 1], Vp is robust and the sequence Vp converges robustly to Vp. 

Proof. By a standard coupling argument, assume that Zk and Zj^ are defined on the same prob- 
ability space rj. Then, letting Y be either Zk or Z]. , for a; G O, let x{t:Y{uj)) be equal to one 
if trajectory Y{uj)^s first jump, starting at time t, is into a state of G at time t' G [t + Ta, i + T5], 
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and zero otherwise. Similarly, let x (t^Y^uj)) be 1 if y(a;)'s first jump, starting at time t, 
is into G^^' at time t' G [t + Ta,t + Th], and zero otherwise. Then P{t) = Klxit, Zk)], and 
P(^)(t) = E[x'^^\t,zlf^)]. It follows that 

|E[x(t,Zfc)]-E[x(^)(t,zf))]| < E[\xit,zk)-X^''\t,zk)\] 



(1) 
+ E[\x^^\t,z,)-x^^Kt,zf^] 

^ V 

(2) 



Consider term (2) above. We can partition trajectories into two measurable subsets: Qi 
{a; G fi I Zk{t,ui) = Z\ {t,uj),t < ti + Tb} and Oq = f^\ J^i. Let fiQ be the probability measure in 



rj. Applying Theorem 3.2 up to time ti + Tf,, we have that x (*> ■^fc {^)) = X {tj ^^{uj)) for 
u £ ill and P(r2o) < £n- Mence, for any t G [to,ti], 






Iflo 
< Eat -> 0. 

Notice that sn does not depend on t. 

Let us focus now on term (1) in the inequality above. Let Ti < T2 < . . . < Tf^he all the points 
in Disc{G) (which are finite in number as G is robust). Fix t G [ioi ^i]- As G^ ' converges robustly 
to G, for N > Nq they differ only in disjoint balls B(Ti,£), for e small enough. Furthermore, if G 
has a discontinuity for state s in Tj, then the value of G on the left of B{Ti,£) is different from the 
value of G on the right of B(Ti, e). 

It follows that the only trajectories of Zk for which xi'^^^k) 7^ X i'^^^k) are those jumping 
within the set Dq (intersected with [t, t+T;,]), ^^ As the rate functions of Zk are piecewise analytic, 
they are bounded by a constant A, thus the probability of a trajectory jumping in Dq is bounded 
by / (AT) Ae'^^dt < f (n) Mt = Afie{Di^^) -^ (independently of t). It fohows that 

^G ^G 

|E[x(t,Zfc)]-E[x(^)(t,zf))]|<5^, 

with 6]\f = en' + ^fJ-iiDQ ) — )• independently of t, which proves uniform convergence of P^ '{t) 
to P{t) . 

Let us turn now to point 2 of the lemma. 
Consider the set Hp of values p G [0, 1] for which either (i) P{t) is constantly equal to p in one 
analytic piece of P, or (ii) P{t) = p and P'{t) = for some t, or (iii) P{t) = p and P is not analytic 



in t. By Prop. 4.1 and the definition of piecewise analytic functions, the set Hp is finite. Fix a 
p Hp. For such a p, the function P{t) — p defines a robust time-varying set, as it has a finite 
number of simple zeros, all in analytic points of P. 

Call A the set of points in which Vp has a discontinuity, which is finite. Fix e and define A;, to be 
Ujg^ i?(t,e), where i?£(t) = {t — e,t + e). Now, if VF is a neighbourhood of A, then for a small 
e > 0, Ae C W. Let /p(t) = |-P(t) — p\ and consider the set Is = I \ A^. Now, I^ is compact 

^^Notice that robustness of G is not necessary for this proof, but we enforce it for uniformity with the convergence 
of reachabiUty probabihties in Section [5j 
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and fp{t) is different from zero in I^, so that min{/p(t) | t £ I^} = nif, > (by the Weierstrass 
theorem [40j). As P^'^> converges uniformly to P, there is Nq such that, for all N > Nq and all 
tele, \P^^\t) - P{t)\ < ^, hence for all iV > iVo and all tele, Vp{t) = vi^\t). It follows that 
Vp {t) converges robustly to Vp{t). ■ 



A. 2 Reachability 



Proposition (5.1). Let X^ > be a sequence of CTMC models, as defined in Section 



Z, and Zk be defined from X^^> as in Section 3.3. Assume that the infinitesimal generator 



3.1 



and let 



matrix Q{t) of z^ is bounded and integrable in every compact interval [0,T]. Then 

Preach{zl ,t,T, G, U) -> Preach{zk,t, T, G, U), uniformly in [io, ii], as TV -^ oo 

i.e. SUPtg[tp_ti] ll-Preac/i(-^fc ' ,t,T ,G ,U) - Preach{zk,t,T,G ,U)\\ -^0. 

he proposition follows from a similar argument to that used in the first part of the proof of 



Proof. 
Lemma 



4.1 



By a standard coupling argument, we can assume that the processes Z\, and z^ are 



defined on the same probability space J7. Therefore, there exists a sequence e^ G IK+, £n -^ 0, such 

that 

f{oj G rj I Vt < r', Zj^ (w,t) = Zk{u:,t)} > 1 — £n- This means that with probability 1 — en, the 

trajectories of the two processes are the same up to time T' . 

Now, we can define a (measurable) function x = Xt,T,G,u on the trajectories of the CTMCs 
which is equal to 1 if they satisfy the reachability property, and otherwise. Therefore, it holds that 
Preachi^k ,t,T,G,U) = K[x(Zj^ )], and similarly for Zj.. With a similar notation as in Lemma 
let Qi = {uj \ Zl '{t,uj) = Zk{t,uj),yt < tQ + T}, VLq = {uj \ z\, '{t,uj) / Zk{t,uj)}, and /Un be 



4.1 



the probability measure in Q, (i.e. in the trajectory space). Observe that xi^k ) ~ xi^k) on $7i 
and P(ilo) < ^N, hence 

|E[x(zf))]-E[x(zfc)]| < EMzi""^) - X(zk)\] 



k 

xiZk ) -x{zk)Wn 



Hi 



\X{ZP) -x{zk)Wn 



+ 

< eTV ^ 0. 

Uniform convergence follows from the fact that the sequence £m does not depend on the initial or 
the final time of the reachability property, if they are both less than T + ti. ■ 



Lemma (5.1). Let X^ ' be a sequence o f CT MC models, as defined in Section 



3.1 



and let Zj^ 

and Zk be defined from X^^' as in Section 3.3. with piecewise analytic rates, in a compact interval 
[0,T'], for T' sufficiently large. 

Let G{t), U{t), t G [to, ti+T\ be compatible and robust time-varying sets, and let G^'^'{t), U^^\t) 
be sequences of time-varying sets converging robustly to G and U, respectively. 
Furthermore, let P{t) = Preachizk,t,T,G,U) and 

P^'^Ht) = Preach{ZP,t,T,G'^^),U^^)), t G [to,t^]. 
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Finally, fix p e [0,1], MG {<,<,>,>}, and let Vp{t) = I{P{t) cxp}, vi^\t) = /{pW(t) cxp}. 
Then 

1. For all but finitely many t G [to,ti], P^ '{t) — t- P{t), with uniform speed (i.e. independently 
oft). 

2. For almost every p G [0, 1], Vp is robust and the sequence Vp converges robustly to Vp. 



Proof. As in the proof of Lemma 4.1| by a standard coupling argument assume that z^ and Zj^ 



are defined on the same probabihty space O. Then, letting Y be either z^ or Z^ , for a; G fi, let 
x{t, Y{u)) be equal to 1 if trajectory Y{uj) satisfies the reachability problem with respect to G and 
U and zero otherwise, starting at time t and x (i^^i^)) be 1 if Y{lo) satisfies the reachability 
problem for G^ ', U^ ', and zero otherwise, starting at time t. Then P{t) = K[x{t, Zk)], and 
P(^)(t)=E[x(^)(*,4^^)],and 

|E[x(t,Zfc)]-E[x(^)(t,zf))]| < E[\xit,Zk)-X^''\t,Zk)\] 



(1) 
+ E[\x^^Ht,z,)-x^^\t,zi:'^] 

^ V 

(2) 



Term (2) above is bounded by e^ -^ 0, as in Lemma 4.1, by a straighforward application of 



Theorem 3.2, Term (1) is also treated similarly to Lemma 4.1, with an extra argument to deal 
with pointwise discontinuities in the reachability probability. Let Ti < T2 < ... < T/^ be all 
the points in Disc{G) U Disc(U) (which is finite as G and U are robust). If we suppose neither 
t nor t + T coincide with one of the previous points (i.e. all discontinuities are internal in the 
time domain), then by robust convergence of G^ ' (resp. U^ ') to G (resp. U), for A^ > Nq they 
differ only in small disjoint balls B{Ti,e) internal to [t,t + T]. Reasoning as in Lemma 4.1 it 
follows that the only trajectories of z^ for which x(t, z^) ^ x {t, Zk) are those jumping within 
the set D^ ' = Dq U D\j Y\ As the rate functions of Zk are piecewise analytic, they are 
bounded by a constant A, thus the probability of a trajectory jumping in D^^> is bounded by 
/^(jv) Ke-^Ht < /^(AT) Adt = A//£(L»(^)) -^ (independently of t). It follows that, if t T^, with 
Td = {Ti, ...,Th,Ti-T,...,Th-T}, then 

|E[x(t,Zfc)]-E[x(^)(t,zf))]|<5^, 

with 6n = eN + A/i^(L»(^)) -^ 0. 

On the contrary, if f G Td, then a discontinuity of G or U happens exactly at the boundary 
of the time domain [t, t + T] in which we have to verify the formula. In this case, the value of 
sets G^ ' and G (or U^ ' and U) may never be the same at this extreme point t*, whatever small 
neighbourhood of t* in [t, t + T] one takes into account (e.g. if t* = Ti is the left extreme of the 
time domain, it may happen that all changes of G^ ' occur before this point). Therefore, there 
can be a set of trajectories of measure > that are accepted by x and refused by x (or vice 
versa). In particular, this can happen if and only if P{t) has a discontinuity in one of those points 



^^If G or U are not robust, then even if they have a finite number of discontinuity points, the previous argument 
may not hold. In fact, they may have a discontinuity point Ti such that Gs{Ti) — 1 but Gs{t) = in a neighbourhood 
W \ {Ti} of Ti. In this case, it is possible that Gl (i) = on all W, which implies that x(i, ^k)) / x'^' (^1 ^k) for 
all those trajectories that are in state s at time T. 
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(otherwise, convergence follows by continuity). Hence, in these time points, convergence may not 
hold. However, the set T^ is finite, hence point 1 of the Lemma is proved. 



Let us turn now to point 2 of the lemma, which is similar to Lemma 4.1 with extra care for 
the discontinuities of P. 



As in Lemma 4.1, construct the set Hp of values p G [0, 1] for which either (i) P{t) is constantly 
equal to p in one analytic piece of P, or (ii) P{t) = p and P'{t) = for some t, or (iii) P{t) = p 
and P is not analytic in t. This set is finite, and for p Hp, the function P{t) —p is easily seen to 
define a robust time- varying set, as it has a finite number of simple zeros, all in analytic points of 
P. Consider now the set A of discontinuity points of Vp. Fix e and define A^ to be UteA-^(^'^)' 



where -Be(t) = {t — e,t + e). By reasoning as in the last part of the proof of Lemma 4.1, letting 
Je = /\^£ and min{|P(t) —p\ \ t G I^} = rris > 0, as P' •*(t) converges in I^ to P{t) with uniform 
speed, there is Nq such that, for all N > Nq and all t £ Is, |P(^)(t) - P{t)\ < ^, hence for all 

iV > iVo and all t G 4, Vp{t) = VJ'^\t). It follows that \4'^\t) converges robustly to Vp{t). 

However, here we need extra care as the set I^ may contain time instants t in which the 
convergence of P^^' to P does not hold, but that do not generate a discontinuity in Vp, because 
P{t~^) and P{t~) are both greater or both less than p. These points do not create problems, 
essentially because the function P^ ' , for A^ large, remains close to P. In fact, convergence at i 
fails because the jumps in C'^-* and U^^^ happen at time instants converging to the ones of jumps 
in G and U, but not necessarily at t. This slightly puts out of synchronization the time at which 
the discontinuity happens, but the values of p(^) and P around such a discontinuity are close. 
This implies that, for A^ large, P'^^' will remain below p if both P{i~^) and P{i^) are below it, 
and similarly for the symmetric case. 

A formalisation of this argument requires a more careful inspection of the behaviour of G^^' 
(respectively U^ ' ) near a discontinuity of G (respectively U) , and a clarification of the connection 
between discontinuities in G and U and discontinuities in P. For the former point, note that by 
the robust convergence property of G' ' to G, if G has a discontinuity for state s at time t, say 
from to 1, then G^^' also has a discontinuity of the same kind near t. In fact, it can do more 
than one jump around t, but for sure, for any small e > and A^ large, it will equal before t — e 
and 1 after t + e. The point is that these additional jumps do not matter, as they happen so close 
to each other that almost no probability mass moves in between, hence they have a vanishing effect 
on p(^) (as N grows). As for the connection between discontinuities in G and U and the function 
P, observe that we can have a discontinuity in P at time t only if either G ox U has a discontinuity 
at time t or at time t + T (they cannot both have such a discontinuity, due to the compatibility 
condition). There are many cases to take into account (a change from goal to non-goal, or from 
non-goal to goal, and so on), but only a few of them induce a discontinuity, specifically a change 
from non-goal to goal of a safe state s at time t + T (inducing a discontinuity in any safe and 
non-goal state at time t), a change from goal to non-goal of a safe state s or form unsafe to safe of 
a non-goal state s at time t (inducing a discontinuity in s), and a change in the goal status of an 
unsafe state at time t. In the first case, we can have a discontinuous increase in P. In the second 
case, the value of P in s can drop from 1 to a value p' < 1. In the third case, the value of P can 
increase from to a value p' > 0. In the fourth case, which is a rather strange case, the value of 
P changes from to 1, or vice versa. 

To understand the connection between P, P^ ' , G and G^ ' , consider a situation of the first kind, 
in which one or more safe states s change from non-goal to goal at time t + T. This creates a 
discontinuity in P{t) for any safe state s' , such that there is a non null-probability of going from 
s' to s along a safe and non-goal path from t to t -|- T. This probability, in fact, is added to 
Pg' (t) , according to the discussion in Section [5] Suppose for simplicity that only a single state s 
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changes status in t + T from non-goal to goal. Then this happens close to t + T also in G^'^K In 
fact, s can change status many times in G^ ' , near t + T, but only the first one really matters 
for the discontinuity of P'^"> . This happens because the probability added to P^, for subsequent 
jumps of s from non-goal to goal state close to t + T is only the probability of jumping into s 
from another safe state in the short time interval in which s is non-goal. To be more concrete, 
if s jumps from non-goal to goal at time t\ + T, then from goal to non-goal at time tg + ^ 
and back to goal at time tg +T, then the probability added to P^, (tg ) is bounded by the 
amount of probability mass that can flow into s in between times t2~\-T and tg + T, which 
is of the order of ig — ^2 ' hence vanishes as N grows (as ^2 ^iiid tg collapse to t). More 
precisely, if A is an upper bound for the exit rate of the single agent (uniform in A^, which can be 
found as the exit rate of Zj, converges to the exit rate of z^, which is itself uniformly bounded), 
then the jump size at time tg is bounded by 2A(t3 — ^2 )• Furthermore, the speed at which 
P^, can increase or decrease, excluding jumps, is also bounded by 2A, so that the value of P^, 
cannot vary too much after the first jump in a small time interval of size A^ around t. In fact, 
combining these two arguments, the total variation (after the first jump) is bounded by 2AAt ^^ 
Note that, if more than one safe state changes goal status at time t + T in G, then in G^^ 
these events can happen asynchronously, hence to see the full increase in P^, we need to wait 
that all those states do their first jump in G^ ' . Yet the bound in terms of A and At remains 
valid. The other discontinuous jump types are treated analogously (with the exception of the 
jump from to 1 or from 1 to 0, which however contains any threshold p in its interior). We 
can now give a formal argument that discontinuities in I^ are not a problem. Assume that P has 
a discontinuity in t such that /x = imix{Ps^{t~^),Ps/(t~)} < p, and call e = p — p.. Now, choose 
5 such that 4(5A < e/4, ||Py(t - 5) - Ps'(*")ll < e/4, and \\Ps'{t + 5) - Ps'{i'^)\\ < e/A. Then, 
choose an Nq such that, for N > Nq, all the jumps of G^^' are closer than 6 to the jumps of 
G, and such that \\P^, {i ± 5) — Ps'{i i: S)\\ < e/4. Then, using the previous reasoning, we can 
see that sup^^]i_s^i^s[Pf\t) < ma^{PiP {i + 6) , P^f\i - 5)} + A6A < p + 3/4e < p. The first 
inequality holds because max{P^, {i + S),P^, {i — 6)} is a value close to the value of p(^) after 

the first jump, and is combined with the bound 46A on the variation. Hence P^, ultimately does 
not cross the line p around t. The case in which uim{Ps'{t~^), Ps'{t~)} > pis dealt with similarly. ■ 

A. 3 CSL model checking 



Theorem (6.1). The CSL model checking for ICTMC, for piecewise analytic interval computable 



rate functions, is decidahle for a robust CSL formula (p{pi, ■ ■ ■ ,Pk)- 

Proof. First of all, we prove that we can approximate the function P{t) for any top next or until 
formula ip with arbitrary small precision. To start, notice that procedures for integrating ODEs 
and doing matrix multiplication (which are at the basis of the methods in Sections H^ and ^ can 
be computed with arbitrary precision, due to the assumptions of interval comput ability. Hence, 
let us focus on the set of zeros of P{t) — p, for a dependent next or until formula (pi. We want 
to prove that we can find those zeros with arbitrary precision, and that in doing this we will be 
able to compute the probability of any next or until formula which contains c/Ji as a sub-formula 
with arbitrary precision. If ip is robust, then the time-varying truth of formula ipi is robust. 



■^°The factor 2 comes from the fact that we are working with a combination of the backward and forward equation, 
both giving an upper bound of A on the rate of change of probabihty mass. 
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This means that P{t) — p has a finite number of simple zeros (i.e. their derivatives are not zero). 
Hence, it is possible to effectively encapsulate them in disjoint intervals of size as small as desired 
p-gpTl Therefore, we can compute the time-varying truth value of the set of states satisfying the 
formula ipi with arbitrary precision, in the sense that for each e > small enough, we can provide 
intervals of size at most e, each containing a single discontinuity point of the set in which one 
or more states change truth status. The condition on compatibility ensures that we can combine 
such approximation of time- varying sets and still obtain robust sets. (This may fail if we take 
the minimum (conjunction) of two truth sets which have a discontinuity for s in the same time 
point T: we can obtain a function which is neither left nor right continuous, a situation that 
cannot originate from a simple zero.) Furthermore we have the further property that we can 



always assume that there is a zero in every approximation interva Consider now the problem 



of computing the probability of an until formula, having two approximations of time- varying truth 



as described above. Reasoning as in the proof of Lemma 5.1 we can see that if we choose an 
arbitrary point in each interval wrapping a discontinuity point in spite of the correct one, we 
commit an error in computing the probability of the until which is uniformly bounded by the 
total size of the approximation intervals. Hence, we can make such error as small as desired. A 
similar conclusion can be drawn for a next formula, invoking the line of reasoning of Lemma |4.1[ 
Reasoning inductively, we can therefore compute with any arbitrary precision the probability -P(O) 
of any top until formula. 

Given this value, we then have to solve the inequality Ps{0) < pi (or -Ps(O) > pi) for any s 
and any top next or until formula ifi. By the robustness of the CSL formula ip, it cannot be 
that Ps{0) = Pi, hence we can effectively solve that problem by computing -Ps(O) with precision 
£i < l-Ps(O) — Pi\- As we are doing interval arithmetic computations, we can increase the precision 
until each pi will be outside the approximation interval for ^,(0). This proves that the algorithm 
presented is effective for robust formulae and eventually computes the exact answer. ■ 



Theorem (6.2). Given a CSL formula (p{p), with p G [0, 1] , then the set {p | f{Tp) is robust} is 



relatively opei] \ in [0, 1] and has Lebesgue measure 1. 

Proof. We will prove the theorem by structural induction on the formula ip. We first need some 
preliminary definitions. Consider an until formula (p = 7^x]p(y'iU^^' ^1(^2) or a next formula 
ip = Ptxjp(X[-'"'"^''](/9i) and call q a generic tuple of values for the thresholds on which ipi and (in the 
until case) ip2 depend on. Fix a q such that the time-varying sets for c^i and (/?2 are robust. An 
open neighbourhood C/q of q is robust if the time-varying sets for ipi and <p2 are robust for each 
q' £ f7q. Observe that in C/q the number of discontinuities of time-varying truth sets of ipi and ip2 
does not change {ipj is robust for each point in [/, and a change in the number of discontinuities 

^'^As the number of zeros is finite and tlieir first-order derivative is non-zero, tlie function /p(t) — P{t) —p crosses 
zero in tiiose points. Furtliermore, notice tliat tiie absolute minimum value of the derivative in those zero points is 
> 0. Hence, there is an e such that each interval of size e containing a zero point has different signs at the extremes 
and the derivative is provably different from zero. By iterated bisection, we can always find such intervals after a 
finite number of steps. Furthermore, all intervals J not containing a zero can be eventually discarded by bisection, 
computing an upper bound L on the absolute value of the derivative in such intervals and bisecting them until we 
can prove that they are disjoint from zero using the Lipschitz condition with Lipschitz constant L (compute fp on 
a single point a; in J of length S, and discard J if |/p(a::)i — LS > 0). 

■^^If we take the minimum (conjunction) of two truth sets which have a discontinuity for s in the same time point 
T, then even if the conjunction is robust, when we have an approximation of the time- varying truth function, we 
can never know if both discontinuities happen in the same time point or in different ones. 

^^A set U G V is relatively open in 1/ C W, where W is a topological space, if it is open in the subspace topology, 
i.e. if there exists an open subset Ui (^ W such that U = V nUi. 
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can happen only at a non-robust point) and the time-instants at which such discontinuities happen 
depend continuously on q. Now, we define the set-valued function 6 : C/q — )• 2''^'^' in the following 
way: Given q', 6(q') is the set of values p which causes the time-varying truth set of ip to be 
non-robust. Therefore, 6(q') contains the values of P.s{t) for which Pg{t) is zero, the values Ps{t~) 
and Ps{t'^) for each non-analytic point t of P, and the values of constant pieces of Pg, plus the 
value -Ps(O). Hence it is finite. 
By possibly restricting C/q, we can also assume that the number of points in 6(q') is bounded by 



|6(q)| in C/c^^ and the value of such points depends continuously on q'. Therefore, the set-valued 

map 6 : C/q —7- 2''^'^' is upper- semicontinuous in C/q, i.e. for each neighbourhood C/fe(x) of 6(x) in 

[0, 1], there is a neighbourhood C/x of x in C/q such that 6(C/x) ^ C/^(x)- 

Given a formula ip = (p{p), p € [0, 1]*^ , we define the set R^:, C [0, l]'^ of all robust thresholds, i.e. 

Po G -R if and only if v'(po) is robust for cp. Hence, our goal is to show that R^ is open and has 

measure 1 for any formula (p. 

We are now ready for the inductive argument. 

Base case: The base case corresponds to (boolean combinations of) atomic formulae, which are 
robust for each p G [0, 1] . 

Boolean combinations: The only non-trivial cases are the conjunction or disjunctions of until or 
next formulae. In these cases, we have to enforce the compatibility condition by guaranteeing 
that the discontinuity times of truth-valued functions are disjoint for each pair of until or 
next formulae. Consider two until or next formulae fi and ip2, and let p = {(li,Pi,(l2,P2), 
where pj is the threshold for formula ipj and qj is the set of constants which fj depends 
on. By inductive hypothesis, the robust sets Rj = Rip. for ipj are open and have measure 1. 
Now, let Pj = Pj{t, qj) be the probability of the until or next path formula in ipj, and fix a 
robust point q = (qi,Pi,q2)- Let g{q) be the set valued function g{q) = P2{{t \ Pi(i,qi) = 
Pi},q2) U 62 (q2), where 62 is the set of non- robust points for ip2, as defined above. The 
set ^(q) contains all thresholds for 932 for which c/?2 is non-robust and all thresholds that 
would make the boolean combination non-robust. By properties of the analytic functions, 
it follows that ^(q) is finite. Hence by arguments similar to the ones above, the function 
g is upper-semicontinuou^2J ^^ ^ neighbourhood U of q. Therefore, letting p2 g{q) and 
^ 1^ 5(q) = a neighbourhood of p2 in [0, 1], we can find a neighbourhood C/ of q such that 
g{U) nV = 0, so that W = U x V is an open neighbourhood of p = (qi,Pi,q2,P2) which 
contains only robust points, which proves that R = R^p is open. 
Furthermore, R is a fortiori measurable. Now, let hji : [0, Ij'^i+'^s _^ |q^ ]^| \^q h^q indicator 



'^^ We have to restrict U to avoid the appearance of further zeros of the derivatives away from current zeros. Note 
also that if a value p £ &(q) corresponds to a non-simple zero at a time io in which the derivative has a maximum 
or a minimum, a small perturbation of q can split it in two, or make it disappear. Just think about raising or 
lowering a curve having a local maximum or minimum with value zero. However, split zeros will be at points ti and 
i2 arbitrarily close to to, and therefore, by continuity of Ps{t,q) with respect to q, the values of Ps(fi,q'), for q' 
close to q will be close to Ps(io, q)- Zeros that disappear are not a problem for semicontinuity, as the empty set is 
contained in any set. Hence, we need to count those points twice in |&(q)|. 

^^The number of solutions of Pi(i,qi) = pi in a sufficiently small neighbourhood Ui of (pi,qi) is constant and 
hence the set- valued function (?i(pi, qi) = {t | Pi(t, qi) = pi} is upper semicontinuous. Furthermore, in a sufficiently 
small neighbourhood U2 of q2, the function P^it, q2) is continuous in q2 for each continuity point t of P2{t, q2). Points 
for which P2(i,q2) is not continuous are covered by b, hence both P2{t'^,<l2) and P2(i~,q2) are in p(q). Now, for 
each neighbourhood V of g{ci), by piecewise analyticity and right/left continuity of P2, we can find a neighbourhood 
U2 of q2 and a neighbourhood Vi of (?i(pi,qi) such that b2{U2) C V and both {p \ p = P2 (t"'' , q2 ) , i G Vi} C V and 
{p \ p = P2 {t~ , q2 ) , i G V\ } C V". Now, by upper-semicontinuity of gi , there is a neighbourhood Ui of (qi , pi ) such 
that gi{Ui) C Vi. It follows that g{Ui x U2) C V, hence g is upper-semicontinuous in q. 
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function of the set R in which the boolean combination (/? of (pi and (p2 is robust. Note that 
R C Ri X R2, and cah R2 the set of thresholds q2 for which sub- formulae of ^2 are robust, 
which is open and has measure 1 in [0, 1]'^^"^ by inductive hypothesis. By Fubini's theorem: 



IJ^iiR) = / hR{qi,pi,q2,P2)fi£{d(ii,dpi,dci2,dp2) 

Jf0,ll'=l+'=2 

/ hR{cii,pi,q2,P2)fJ'iidp2)fJ-i{dcii,dpi,dq2) 
J[o,i] 

^ He{dqi,dpi,dq2) = / fJ.e{dcii,dpi) / ^ fie{dq,2) = 1, 



RixR'2 

which proves that R has measure 1. 

If we have a boolean combination of j > 2 until formulae, we simply reason pairwise and 

then take the intersection of the so-obtained robust sets, thus getting an open set of measure 

1. 

Until formulae: Let ip = Vt,<3p^{(pi\J'^^''^'^'ip2)- By inductive hypothesis, the set R^-^ x R^^ ^ 
[0, l]'^"^ for which ipi and 992 are robust is open and has measure 1. By reasoning as in the 
boolean combination case (and considering all until and next conjunct/disjuncts of c^i and 
1P2), we can immediately conclude that the set R' Q R^^ x R^p^ in which the time varying 
sets of 931 and ip2 are robust and compatible is open and has measure 1. 

Now, fix a point q € i?', let ?7 C fi' be a robust neighbourhood of q, and consider the set 
valued function 6 : C/ — t- 2''^'^' as defined above. Now fix p i>{^)l and choose a neighbourhood 
V oi p such that V n 6(q) = 0. As 5 is upper-semicontinuous, there exists W <ZU such that 
hiW) n y = 0, hence if is robust in T^ x y. By the arbitrary choice of p = (q,p), it follows 
that R = R^p is open, and hence measurable. Now, let h^ : [0, l]'^ — ;■ {0, 1} be the indicator 
function of the set R in which ip is robust. By Fubini's theorem, it follows that R has measure 
1. 

Next formulae: The argument for a next formula ip> = Vt^p,^{Xy"-''^''^(pi) is essentially the same 
as for until formulae, with the only difference that the inductive hypothesis is applied only 
to ipi and there is no need to ensure compatibility. ■ 



Theorem (6.3). Let X^^> he a sequence ofCTMC models, as defined in Section 



and Zj^ be defined from X^'^' as in Section 3.3 



3.1 



and let Zj^ 



Assume that Zj^ , z^ have piecewise analytic infinitesimal generator matrices. 

Let (p{pi, . . ■ ,Pk) be a robust CSL formula. Then, there exists an Nq such that, for N > Nq and 

each s £ S 

s,OI= (iv) 99 <^ s,0^zk V>- 

k 

Proof. We use structural induction to prove that, for each formula 99, the time- varying truth sets 
Vip of 99 in Zj^ converge robustly to the robust time- varying truth set V^ of (p in z^- 

Base case: The case for atomic propositions is trivial, as V^ and Vip are constant and equal. 



Negation: Let ip = ->ipi. The result follows because V^^(t) = 1 — V^-^ and V^ (t) = 1 — V^ 



VI 
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Conjunction/Disjunction: Let ip = (pio 932, ° G {A, V}. Due to the compatibility condition of 
robustness of </? with respect to z^, the set V^p{t) = mm{V^p^{t) ^Vip^it)} ^ mm £ {min, max} 
is robust, with Disc{V^) C Disc{Vy,-^) U DiscCVip^). Using the inductive hypothesis, it easily 
follows that V^ (t) = m,m,{Vip-^ (*))H'2 (*)} converges robustly to V^{t). 



Next: Let ip = Vt^p(X\ "" '''(pi). By inductive hypothesis, we can apply Lemma 4.1 and deduce 



that Vip is robust and Vip converges robustly to V^. 
Until: Let 99 = 'P(^.p{ipi\J'''^°"'^''^ (P2) ■ By inductive hypothesis (and the compatibility condition in 



the definition of robustness of f) , we can apply Lemma 
and Vtp converges robustly to V^. 



5.1 



16 



and deduce that V^ is robust 



The fact that Vip converges robustly to the robust set Vtp, combined with property 1 of robust- 
ness of if, let us conclude that the truth value of if at level N converges to the truth value of 
the limit ICTMC at time zero (if was a point in which convergence of probability fails, then a 
small perturbation in p could change the truth value of 99 in the limit ICTMC, contradicting the 
robustness of ip; furthermore, robustness of ip forbids that Ps(0) = p)- ■ 



A. 3.1 Comparison of CSL model checking for Z^ and (Z^, ,X^ ) 



We will prove now the lemmas in Section 6.3 of the p ape r. W e wi ll start by an auxiliary result, 



which is needed to adapt the proof style of Lemmas 4.1 and 5.1 to the processes {Zj^ '-^L ) 



and (zjt,Xfc) discussed in this section. In particular, in Lemmas 4.1 and 5.1 we used the fact that 
processes jump with a small probability in a small temporal neighbourhood of the discontinuity 
points of time-varying sets. In the space-based setting, however, we need to consider neighbour- 
hoods of the boundaries of goals and unsafe sets. Therefore, to use the same proof style, we need 
to bound the time trajectories spend in such a neighbourhood, in a uniform way in space. The 
key point is that, as the convergence results we are interested in depend only on a neighbourhood 
of the trajectory <I>([0,T],xo), we can always choose a small flow tube such that the velocity with 
which a trajectory crosses the boundary of a goal or an unsafe set is close to that of '&([0, T], xq), 
and so will be the time spent in a neighbourhood around such boundary. In the following, we will 
make this intuition formal. 

First of all, observe that the notion of robust set V implies that when a trajectory crosses the 
boundary dV in a point x, the function h defining the smooth manifold of the d-set dV around x 
changes sign. 

We will now prove an upper bound for the time spent by a trajectory in a neighbourhood of the 
d-set D = dG of a robust set G in 5 x £'0, a T, eo-Aow tube of <I>([0, T], xq). 

For each trajectory (^{\r~ {Eq) ,T^ {Eq)],:^.) , x G £^01 consider the points Disc{s,x.) = {(s,xi) G 
Z) I xi G <^{[T-{Eo),T+{Eo)],x)} in which it intersects D. 

We define the e-neighbourhood of D m S x Eq as D^ = \J(sx)eDn(SxEo) ^^i'^'^)^ which is an 
open set. Note that Ds = U{s,x)e5x£'o U(s,xi)eDisc(s,x) -Se('SiXi), as S x Eq is the union of a set of 
trajectories. 



■^''We need to apply it twice for the two reachability problems involved in computing the probability of an until 
formula, noticing that the probability of the path formula within ifi is an analytic combination of the two so-computed 
probabilities. Robustness of Vip and robust convergence of V^ to Vip follows from the same arguments of Lemma 



5.1 



Alternatively, one can modify Lemma 5.1 and tailor it to the reachability involved in the until case (which 
reduces the time window in which one can reach the goal set), by a straightforward modification of the definition of 
x'^' and X S'lid adaptation of the arguments for robust convergence. 
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Now, by the robustness property of G in 5 x Eq, we have that \Disc{s,x.)\ < k. Furthermore, 
by the robustness of G, the trajectory $([0,T],xo) will cross D moving from the interior of G to 
the interior of its complement, or vice versa, for any point (s,x|) S Disc{s,:>CQ) and any s € S. 
Consider a neighbourhood W of (s,x|) in which D = dG is a smooth manifold. Therefore, there 
is a sufficiently smooth function h va.W such that D r\W \s the zero set of h. By the robustness 
property of G, the function /i($(t,xo)) equals in U (the time such that $(t,xo) = x|), and 
changes sign around tj (i.e. it is positive in tj — 5 and negative in tj + 5, for a 5 > 0). It follows that 
the derivative of /i($(t, xq)) in ti is non-null. As it equals V/i(x|) •F(xf ), we have that Vh{x.f) / 0, 

and hence I iig^ggl, •F(xf)| > 0. 

Now, by choosing a suitably small neighbourhood Wi <ZW oi {s^'H-D (we need to ensure that the 
manifold containing (s,x|) is the closest one in Wi, among those constituting D), we obtain that 
the function p{t,x) = dist{{s,^{t,x)),dG) = mf(^sy^^g(j\\^{t,x) — y|| is differentiable (in t and 

x), and its derivative in (0,xf) is /9'(0,x|) = g^^ = iw/hCx")]} ' -^(-^f)' ^-^^ ^^ ^^ ^^^ projection 
of the vector field along the normal to the surface {h = 0}. Now, by continuity of p' , we find a 
neighbourhood W2 C Wi of {s,xf) such that |/o'(0,x)| > /of/2, where p^ = \p'{0,xf)\. 

Now, choose ei < Eq and e that satisfies: (i) the number of intersections between D and a 
trajectory in 5 x £^1 is constant and equal to the number k of intersections of 5 x <I>([0, T], xq) with 
D (this is possible because the flow tube 5 x £^1 is a small neighbourhood of 5 x $([0, T], xq)), and 
(ii) Dg, the e-neighbourhood of D in the T, ei-flow tube S x Ei, is contained in the neighbourhood 
W2 of (s, x|) identified above, for any s £ S and i < k. 

By the choice of W2, it follows that the speed at which each trajectory of S x Ei travels in Dg 
(with respect to the distance from D) is bounded below by po = mins^i pf/2, and hence the total 
time Tg a trajectory of 5 x Si spends in Dg is bounded above by ^, as it has to travel a total 
distance of 2ek with speed no less than pQ. Notice that this bound is independent of the specific 
trajectory considered. 

With the previous discussion, we have proved the following 

Lemma A.l. Let Eq C E be a T,eQ-flow tube for xq. Let G be a robust subset of S x Eq. Then, 
there are positive constants £1, e, and pQ such that, for any e' < e, the total time r^/ a trajectory 
in S X El (El the T,ei-flow tube for xq) spends in Df,i , the e' neighbourhood of D = dG, satisfies 
T^i < -^, where k is the number of intersections of any trajectory with D. 

Equipped with this lemma, we can now prove the following one. 



Lemma (6.1). Let Eq G E be a T,eQ-flow tube for xq. Let G be a robust subset of S x Eq, and 
G^ ' be a sequence of subsets of S x Eq that converge robustly to G. 

Let P{s, x) = Pnextiy, s, X, Tq, Tfo, G) be the probability that the first jump of y{t) is into a state 
in G and happens at a time t £ [Ti, T2], given that y started at time t = in state (s, x) G 5 x Eq, 
and let P^^' {s,x) = P^^J^(Y''^',s,u^^'{x),Ta,Th,G^^') be defined similarly, withG andx replaced 
by G^ ' and u^ '(x), respectively. 
Furthermore, define V = {(s,x) | P(s,x) \xi p} and 

y(^) = {(s,x) I P^^\s,x) txip}. Then there exists ei > such that, in Ei, the (T — Th),£i-flow 
tube for xq: 

1. P' •*(s,x) — )• P{s,x) for all x G Ei, uniformly in (s,x). 

-2. IfVx_g{t), t £ [r~j(E'i),T3+j(£'i)], is a robust time-varying set, then V is robust in Ei and 
V^ ' converges robustly to V. 
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Proof. First, notice that we can always restrict to an arbitrary small neighbourhood of <l>(t,xo), 
i.e. to a T, ei-flow tube Ei for xq, with £i as small as desired. This follows from the convergence 
in probability implied by Kurtz Theorem |3.1[ Given 6 > 0, this guarantees that we can find an 
index A'^o such that, for any N > A'^o, with probability at least 1 — 6, the trajectories of X^ -'(t) are 
contained in Ei. Furthermore, we can choose such an index Nq independently of x. Hence, given 
£"0, if we consider a flow tube Ei for xq with radius ei < eo/2, each flow tube of radius ei wrapping 
a trajectory in Ei will be contained in Eq. This guarantees that the next-step probability for Y' •* 
for any point in S x Ei will ultimately depend only on the goal sets G^^^ within S x Eq. 

We first prove point 1 of the lemma in a (T — Tfe),ei-flow tube Ei for xq, for an £1 < £o/2 to 
be fixed in the following. We will prove convergence of P'^) to P for each (s,x) £ S x Ei. 



We will now use an argument similar to the one of Lemma 4.1 Couple y and Y^^-* on the same 
probability space Q,, and let x (resp. x ) be random variables defined on sample trajectories and 
equal to one if the trajectory's first jump is in G (resp. G^^'). Then, as P(s,x) = E[x(s,x, y(cj))], 
where y(0) = (s,x), and similarly for P(^-*(s,x), to show convergence we just need to prove that 
|E[x(s,x,y)]-E[x(^ns>x,YW)]| ^0. It holds that: 



|E[x(s,x,y)]-E[x(^)(s,x,YW)]| < E[|x(s,x,y) - x^^^l 



s,x,yj 



(1) 
+ E[|xW(.,x,y)-x(^ns,x,YW) 



(2) 



To treat term (1), invoke Lemma A.l assume ei is smaller than the one required by the lemma, 



and let e and po the other two constants obtained from it. Now, as in Lemma 4.1, observe that 
for each e' < e, the only trajectories of y for which x and x can have a different value are 
those jumping at a time at which y(i) is in D^i. Now, the total amount of time y spends in D^i is 



uniformly bounded by r^/ < ^— , where k is the number of intersections of a trajectory in S x Ei 



po 



with D. It follows that term (1) can be bounded by -^ — , where A is an upper bound for the 
jump rate of z^ in S x Eq. 

The bound on term (2), instead, follows from the convergence of Y' •* to y, but it requires 



a slightly different treatment than in Lemma 4.1, as now the time varying sets for Y*-^-* depend 



on the sample trajectories of X' •*, hence they are random quantities. Call G ^(jv)(0 the time- 
varying sets relative to G^^', but defined with respect to trajectories of Y*-^)(t). The time varying 
sets for y, with respect to G^^', are denoted by Gx (t), while that relative to G is Gx(t). We 
will need now to control two things: first, we will construct a neighbourhood of D in such a way 
that all the time varying sets are the same outside it, for A^ large enough. Then, we will bound 
the time taken by X*^ •'(t) to cross such a neighbourhood (again for A^ large enough). 

Assume e' < e/2, and consider the e'-neighbourhood -D^/ of D. Invoking robust convergence, 
choose A^o such that for A^ > Nq, G^ ' coincides with G outside D^/. We now want to find a 
neighbourhood [t — T',i + t'] of the time t in which x(t) G D, such that we are guaranteed that 
if t falls outside this neighbourhood, both x(t) and X^ ■'(t) are outside D^i. For any such time t, 
it clearly holds that G ^.j^Jt) and Gx (t) coincide. To find such neighbourhood of i, let A^i be 



X 



,XW 



such that, for A^ > A''i, ||X(^)(t) — x(t)|| is less than e' with probability 1 — 6 {6 tohe fixed later). 
Call Q,£' this event. Condition on it and consider -D2e'- If x(i) -D2£'; then it follows that X'^(t) 
will not belo ng to D^i. Hence, we just need to bound the time T2£' that x(t) spends in D2e'- By 



Lemma 



A.l 



this time is no more than ^^. 

po 
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Now, using Theorem 3.2 choose an A'^2 such that, for A^ > A'^2j Zj^ (t) coincides in [0, T] with 
Zfc(t) with probabihty at least 1 — 6. To bound term 2, observe that if Zj^ and z^ are the same, 
and conditional on event O^/, if both Zj, and z^ jump at time instants in which x(t) is outside 
D2e') then x {s,'^,y) and x (s,x, Y^^') will have the same value. 

Therefore, we can bound term (2) by the probability of x {s,'^,y) 7^ x' '(s,x, Y*^^), which 
is itself bounded by 

F{zk jumps in L»2.'} + ^{^''1 + H^f^ + ^fc} < ^^ + 2<5. 

Now, fix e > and choose e' < min{e/4, j^^}, and 6 < e/4. By combining the bounds on term 
(1) and term (2), we obtain that 

limsup|E[x(5,x,y)] -E[x(^)(s,x,yW)]| < ^^ + 26 <' + ' =e, 

which by the arbitrariness of e implies that 

hm |E[x(s,x,y)]-E[x(^)(s,x,Y(^))]| =0. 

Therefore, we obtain that P*-^-*(s,x) — )• P(s,x), and this convergence is uniform with respect to 
(s,x) G 5 X ii^i, as the bound derived above is independent of it. 

As for point 2 of the lemma, observe that by the fact that the time- varying set l^Q(t) as- 
sociated with the fluid trajectory <l>(t,xo) is robust, and by piecewise analyticity of P, we can 
choose an £\ sufficiently small not only to satisfy the constraints to derive convergence discussed 
above, but also such that all trajectories in the flow tube E\ are robust, i.e. their time varying 
set with respect to P are robust (just observe that the function P(s, <l>(t,x)) is piecewise analytic 
in t and x for each s). Furthermore, we have chosen e\ so that the number of intersections of 
$(t, x) with V in each state s, i.e. the number of times -P(s, <I>(t, x)) — p changes sign, is the same 
as that of <l>(t, xq). Now, consider the boundary dV in 5 x £'1, which is the zero set of the function 
/i(s, x) = P{s,x) —p. By continuity of P and by the robustness property of Fxo (0 > ^^ have that the 
trajectory (s, <I>([0,r],xo)) intersects dV in points x| in which the function P is analytic. Hence, 
P will be analytic in a neighbourhood of xf , and, by a suitable choice of ei, P will be analytic in 
the whole component of dV containing x.f. It follows that dV is the union of smooth manifolds 
(analytic in this case). 

It follows that, by choosing £1 suitably small, dV is a d-set and V is robust. 

As for the robust convergence of V^ ' to V, by the uniform convergence of p(") to P outside an 
open neighbourhood of dV, we obtain the robust convergence of V^^' to V. ■ 



Lemma (6.2). Let Eq C E be a T,eo-flow tube for xg. Let U and G two robust and compatible 
subsets of S X Eq, and U^'^', G^'^' be sequences of subsets of S x Eq that converge robustly to U 
and G, respectively. 

Let P{s, x) = Preachiy-, s, X, T^, Tfc, U, G) be the probability that y (t) reaches a state in G within 
time [Ta,Th], avoiding any unsafe state in U, given that y started at time t = in state (s,x) G 
S X Eq, and let pW(s,x) = PJ:^J^,^{Y(^\s,u'^^\x),Ta,Tb,U^^\G^^^) be defined similarly, with 
G, U, X replaced by G^ ' , [/'•*, and i^' ■'(x), respectively. Furthermore, define V = {(s,x) | 
P(s,x) M p} and V^ ' = {(s,x) | P^ ■'(s,x) ixip}. Then there exists ei > such that, in Ei, the 
(T — Tij), El-flow tube for xq: 
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1. P'^)(s,x) —7- P{s,'x) for all Ji. £ Ei\ D, where D is a d-set, uniformly in (s,x). 

2- IfV-x-oit), t £ [T^^^{Ei),T^^{Ei)], is a robust time-varying set, then V is robust in Ei and 
V^ ' converges robustly to V. 



Proof. First, notice that, as in Lemma 6.1, we can always restrict on an arbitrary small neigh- 
bourhood of ^{t, xq), i.e. on a T, ei-flow tube Ei for xq, with ei as small as desired, implying that 
the reachability problem for Y^^ for any point in S x Ei will eventually depend only on the goal 
sets G(^) and C/(^) within S x Eq. 

We first prove point 1 of the lemma in a (T— T^), ei-flow tube Ei for xq, for an ei < eo/2 (ei will 
be fixed in the following). Consider the set D in Eq, D = Disc{G)U Disc{U)U^~'^{Ta, Disc{G)) Li 
^'^{Ta,Disc{U)) U ^-^{Tb,Disc{G)) U ^~^{Tb,Disc{U)), containing the discontinuity points of 
G and U and all points that are mapped by the flow to Disc{G) U Disc{U) after Ta or Tf, units 
of time. D is easily seen to be a d-set. In fact, it is closed and it intersects each trajectory a 
finite number of times, as Disc{G) and Disc{U) are d-sets. Furthermore, each smooth manifold 
of G or U, defined as the zero set of the function /i(x), will be mapped by <l>~^(Tj, •), j = a,b, 
into the smooth manifold defined by the function /i($(tj,x)). This function is smooth as <l>(t,x) 
is piecewise analytic and it is at least of class C^. 



Differently from Lemma 6.1, we will prove convergence of P'^"' to P for each (s,x) G (5 x 
Ei)\D. 
Consider now the set Dq = Disc{G) U Disc(U), and define the e- neighbourhood D^ of it as done 



in Lemma 6.1 It clearly holds that D^ — )• Dq, as e — )• 0. 

We will now use an argument similar to the one of Lemma 5.1 Let x (resp. x ) be random 



variables defined on sample trajectories and equal to one if the trajectory satisfies the reachability 
problem of the Lemma with respect to G, U (resp. G^^', U^^'). Then, as P(s,x) = E[x(s,x, y(w))], 
where y(0) = (s,x), and similarly for P^ '{s,'x.), to show convergence we just need to prove that 
|E[x(s,x,y)]-E[x(^Hs,x,Y(^))]| ^0. It holds that: 

|E[x(s,x,y)]-E[x(~)(s,x,YW)]| < E[|x(.s,x,y) - x^^^H^-x.y)!] 

(1) 
+ E[|xW(.,x,y)-x(^)(5,x,YW)|], 



(2) 



From Lemma |A.1[ we obtain constants ei and e that bound the size of the flow tube Ei, and 
of the Df/ n eighbourhood of Dq. Under these constraints, we can reason exactly as in the proof of 



Lemma 



6.1 



to bound terms (1) and (2) by — \-26, where 6 and e' < e/2 can be chosen arbitrary 

small for A^ large enough, concluding that |p( •'(s,x) — P(s,x)| converges to zero, uniformly in 
{S X El) \ D. 

As for point 2 of the lemma, robustness of V follows by the same argument as Lemma |6.1[ 
Notice that Disc{V) is closed, as it is the union of the zero sets of continuous functions (the ana- 
lytic pieces of P), plus the subset Dp of discontinuity points of P such that liminf P(s, x) < p and 
limsupP(s,x) > p, which is also closed. Furthermore, by robustness of Vugi^t), we can choose ei 
such that all points in Dp satisfy liminf P(s,x) < p and limsupP(s,x) > p (we need this because 
Vx,{t) has to be robust for all x in 5 x Ei). For the robust convergence of V^'^^ to V, Disc{V) 
is a d-set, hence we can use uniform convergence of P^^' to P outside an open neighbourhood 



of Disc{V). Additionally, notice that, as in the proof of Lemma 5.1, the points in which we do 



not have convergence of P^^> to P and that are not in Disc{V), do not create problems, as in a 
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small neighbourhood of those points, P is always strictly above or below p, and the lim sup or the 
lim inf of P^^' in those points will uniformly satisfy the inequality defining V^ ' (thanks to the 
compatibility condition of G and U) . ■ 



Lemma (6.3). Let X^ > be a sequence of CTM C models, as defined in Section 
and Zj^ be defined from X^'^' as in Section 



3.3 



3.1 



and let Zj^ 



Assume that there is a flow tube Eq of xq such that all trajectories in Eq are piecewise analytic. 
Let 9? = V3(p) be a robust CSL formula for the trajectory <l>(t,xo). Then, there is an Nq such that, 
for all N>No, 

S,Xo \=y (f <^ S,l^^^\xo) \=Y(N) (p. 

Proof. We will prove by structural induction on the formula ip, that there is a T, e-flow tube E^p 
of Xq such that V^ converges to V^ robustly, where V^ is the set of points (s, x) G 5 x E^ such 
that s, X |=y If and V^ is the set of points (s, x) G 5 x E^ such that s, i^^^-'(xo) |=y{jv) ^■ 

Base case: the result for atomic formulae is trivial as their truth value depends only on s, hence 
we can choose E^ = Eq and V^ = V^ = S^p x Eq, where 5^ = {s | s |= ip}. 

Negation: If (/9 = ^(fi, we can choose E^, = E^p^, and simply observe that robust convergence of 
Vp^ to F^j implies robust convergence of V^ = E^^ \ V^^ to 1^^ = E^^ \ V^-^ . 

Conjunction and Disjunction: If ip = ipi o p>2, o g {A, V}, consider E^^, a T, Sj-flow tube, and 
sets F^. — )• K/,. . As p is robust for the trajectory starting in xq and xq G E^., by the 
compatibility condition of (/? there exists an e such that the d-sets of V^-^ and V^j ^^^ disjoint 
in 5 X E^p, for E^p the T, e-flow tube in xq. It easily follows that V^ = V^-^ • V^^ is robust in 

5 X ^^, • G {n, U}, and V^P • FJf ^ -^ V^, • V^^ robustly. 

Next: If ip = Vt>^p(X.'-'^'^''^'^'ipi), let E^-^ be a T, e-flow tube for ipi and let V^-^ be a robust set, 
such that Vp-^ — )• V^p^ robustly. By considering the e,r-flow tube Ej for xq, and using the 



robustness of ip, we satisfy the hypothesis of Lemma 6.1, hence there is an (T — T2),e-fiow 



tube Ep for xq such that V^ is robust in S x E^ and V^ ^Vp robustly. 
Until: If (p = V^p{(pi\J^^^''^'^^ip2), let Ep,^ be T, Sj-flow tubes for tpi, i = 1,2, and robust sets 



Vp^, such that V^^ — )■ V^. robustly. By letting £q < min{ei,e2}, such that dVp^ n dVt^ 



'f2 



(which can be found by the compatibility condition enforced by robustness of ip), considering 
the eoi^-flow tube Ej for xq, and using the robustness of ip, we satisfy the hypothesis of 
Lemma 6.2 , hence there is an (T — r2), e-flow tube E^p for xq such that V^p is robust in 5 x E^p 



and Vp — )• V^ robustly. 

Given a formula ip, and the flow-tube E^p for xq, such that V^ is robust in 5 x E^, and Vp — >• V^ 
robustly, then the lemma follows by observing that, due to robustness of tp, xq does not belong to 
the d-set dV^p, hence there is an Nq such that, for all A^ > Nq, (s,xo) G V^ <^ {s,'>^o) G V^p. ■ 
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